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ABSTRACT 

The  YSOVAR  (Young  Stellar  Object  VARiability)  Spitzer  Space  Telescope  observing  program  obtained  the  first 
extensive  mid-infrared  (3.6  and  4.5  pm)  time  series  photometry  of  the  Orion  Nebula  Cluster  plus  smaller  footprints 
in  1 1  other  star-forming  cores  (AFGL  490,  NGC  1333,  Mon  R2,  GGD  12-15,  NGC  2264,  LI 688,  Serpens  Main, 
Serpens  South,  IRAS  20050+2720,  IC  1396A,  and  Ceph  C).  There  are  ~29,000  unique  objects  with  light  curves  in 
either  or  both  IRAC  channels  in  the  YSOVAR  data  set.  We  present  the  data  collection  and  reduction  for  the  Spitzer 
and  ancillary  data,  and  define  the  “standard  sample”  on  which  we  calculate  statistics,  consisting  of  fast  cadence 
data,  with  epochs  roughly  twice  per  day  for  ~40  days.  We  also  define  a  “standard  sample  of  members”  consisting 
of  all  the  IR-selected  members  and  X-ray-selected  members.  We  characterize  the  standard  sample  in  terms  of 
other  properties,  such  as  spectral  energy  distribution  shape.  We  use  three  mechanisms  to  identify  variables  in  the 
fast  cadence  data — the  Stetson  index,  a  x2  fit  to  a  flat  light  curve,  and  significant  periodicity.  We  also  identified 
variables  on  the  longest  timescales  possible  of  six  to  seven  years  by  comparing  measurements  taken  early  in  the 
Spitzer  mission  with  the  mean  from  our  YSOVAR  campaign.  The  fraction  of  members  in  each  cluster  that  are 
variable  on  these  longest  timescales  is  a  function  of  the  ratio  of  Class  I/total  members  in  each  cluster,  such  that 
clusters  with  a  higher  fraction  of  Class  I  objects  also  have  a  higher  fraction  of  long-term  variables.  For  objects  with 
a  YSOVAR-determined  period  and  a  [3.6]— [8]  color,  we  find  that  a  star  with  a  longer  period  is  more  likely  than 
those  with  shorter  periods  to  have  an  IR  excess.  We  do  not  find  any  evidence  for  variability  that  causes  [3.6]— [4.5] 
excesses  to  appear  or  vanish  within  our  data  set;  out  of  members  and  field  objects  combined,  at  most  0.02%  may 
have  transient  IR  excesses. 

Key  words:  circumstellar  matter  -  stars:  pre-main  sequence  -  stars:  protostars  -  stars:  variables:  general 
Online-only  material:  color  figures 
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1.  INTRODUCTION 

Optical  variability  was  one  of  the  original  defining  character¬ 
istics  of  the  class  of  object  later  determined  to  be  stars  in  the 
process  of  formation  (Joy  1945;  Herbig  1952),  or  young  stel¬ 
lar  objects  (YSOs).  Optical  and  near-infrared  (NIR)  monitoring 
over  timescales  of  weeks  to  months  of  the  nearest  star-forming 
regions  (SFRs)  have  shown  that  the  surfaces  of  YSOs  are  often 
mottled,  with  both  hot  spots  (where  gas  accretion  columns  from 
the  inner  disk  impact  the  stellar  surface)  and  cool  spots  (starspots 
analogous  to  sunspots;  Rydgren  &  Vrba  1 983;  Vrba  et  al.  1986; 
Bouvier  et  al.  1993).  Because  the  stars  are  also  rotating,  the 
presence  of  spots  causes  their  apparent  luminosities  and  colors 
to  vary  with  the  stellar  rotation  period.  As  summarized  in  Herbst 
et  al.  (1994),  cool  spots  are  found  on  YSOs  without  disks,  or 
at  least  without  substantial  accretion  disks  (weak-lined  T  Tauri 
stars,  or  WTTs),  which  is  expected  since  those  stars  do  not  gen¬ 
erally  have  other  signatures  of  active  accretion;  however,  both 
cool  and  hot  spots  have  been  identified  on  YSOs  with  substantial 
disks  (classical  T  Tauri  stars,  or  CTTs).  The  largest  amplitude, 
most  variable  optical  light  curves  are  generally  attributed  to  hot 
spots  (Vrba  et  al.  1993). 

The  periodicities  found  in  spot-dominated  light  curves  have 
been  taken  to  be  the  rotation  period  of  the  star,  and  the  derivation 
of  periods  has  long  been  the  most  common  analysis  of  time 
series  data  of  young  stars.  For  solar  mass  YSOs  and  ages  ~ 
few  Myr,  the  distribution  of  rotational  velocities  is  bimodal, 
with  one  set  of  stars  having  periods  on  the  order  of  2^1  days  and 
the  other  with  characteristic  periods  of  8-12  days  (e.g.,  Cieza  & 
Baliber  2007,  and  references  therein).  This  period  distribution 
has  been  interpreted  in  terms  of  a  model  where  the  rotational 
periods  of  the  accreting  stars  are  magnetically  locked  to  the 
Keplerian  rotation  period  of  their  inner  disks  (with  periods  on 
the  order  of  10  days),  whereas  stars  that  are  no  longer  accreting 
spin  up  as  they  contract,  thus  associating  the  short-period  peak 
in  the  rotation  period  distribution  with  stars  that  have  lost  their 
disks  at  young  ages  (e.g.,  Bouvier  et  al.  1997).  This  correlation 
at  young  ages  appears  to  persist  to  later  ages,  with  the  slow 
rotators  on  the  zero-age  main  sequence  (the  stars  with  long- 
lived  accretion  disks)  being  more  likely  to  have  debris  disks, 
which  could  suggest  that  these  slow  rotators  are  more  likely  to 
have  formed  planets  (see,  e.g.,  Bouvier  2008;  McQuillan  et  al. 
2013a,  2013b). 

In  the  past,  periodicities  have  been  most  frequently  deter¬ 
mined  using  ground-based  optical  time  series  observations.  Op¬ 
tical  observations  are  primarily  sensitive  to  phenomena  asso¬ 
ciated  with  the  stellar  photosphere  or  with  other  energetically 
“hot”  regions  (hot  spots,  accretion  columns,  chromospheres), 
and  are  limited  in  regions  of  high  extinction.  In  contrast,  ob¬ 
servations  at  longer  wavelengths  penetrate  extinction  and  also 
offer  a  new  perspective  by  being  sensitive  to  variability  associ¬ 
ated  with  “warm”  or  “cool”  regions — the  disks  and  envelopes 
of  YSOs.  The  dominant  contributions  to  YSO  photometric  vari¬ 
ability  in  the  IR  include  the  Rayleigh-Jeans  tail  of  the  hotter 
processes,  as  well  as  dust  reprocessing  of  emission  from  these 
hotter  processes,  along  with  phenomena  uniquely  associated 
with  the  disk.  Relevant  disk  processes  might  involve  thermal 
emission  from  an  overdense  (or  overwarmed)  region  of  the  in¬ 
ner  disk,  variable  disk  accretion,  structure  in  the  disk  rotating 
into  and  out  of  view  causing  changes  in  the  measured  Ay  toward 
the  star,  or  disk  instabilities  (e.g.,  Fedele  et  al.  2007;  Plavchan 
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et  al.  2008a,  2013;  Herbst  et  al.  2010).  Finally,  standard  geomet¬ 
ric  effects  due  to  orbiting  companions  can  also  be  probed  in  the 
mid-infrared,  uniquely  so  for  more  embedded  sources.  Because 
many  more  physical  processes  can  affect  the  variability  of  YSOs 
in  the  infrared,  relatively  few  infrared  light  curves  are  periodic 
and  thus  straightforward  to  analyze  (see,  e.g.,  Cody  et  al.  2014; 
Morales-Calderon  et  al.  2011). 

One  of  the  first  monitoring  programs  of  YSOs  at  NIR 
wavelengths  was  Skrutskie  et  al.  (1996),  which  monitored 
15  YSOs  in  Taurus-Auriga.  They  found  periodic  variability, 
as  well  as  variability  due  to  accretion  and  extinction. 

The  first  large  program  of  time  series  photometry  of  YSOs 
at  wavelengths  longward  of  1  micron  was  by  Carpenter  et  al. 
(2001,  hereafter  CHS01).  CHS01  obtained  JHKS  monitoring  of 
~3  deg2  of  the  Orion  Nebula  Cluster  (ONC)  over  a  ~1  month 
time  period  as  part  of  the  Two-Micron  All  Sky  Survey  (2MASS; 
Skrutskie  et  al.  2006).  About  1000  Orion  members  showed 
NIR  photometric  variability  in  their  data.28  Typical  light  curve 
amplitudes  were  of  the  order  of  0.2  mag,  but  with  some  stars 
exhibiting  amplitudes  up  to  2  mag;  periods  were  determined  for 
about  a  quarter  of  their  stars.  CHS01  attributed  the  variability 
for  somewhat  more  than  half  of  the  stars  to  cool  spots;  they 
suspected  that  most  of  the  others  could  be  explained  by 
hot  spots,  variable  extinction,  or  variable  accretion.  However, 
they  could  not  make  a  definitive  determination  and  suggested 
multiple  mechanisms  could  be  involved.  2MASS  monitored 
other  SFRs  as  well,  including  Chamaeleon  and  Rho  Oph. 
Carpenter  et  al.  (2002)  reported  on  the  more  limited  2MASS 
study  of  the  Chamaeleon  SFR  and  similarly  characterized 
variability  amplitudes  and  behaviors,  along  with  identifying 
new  candidate  young  star  members  via  their  infrared  variability. 
Plavchan  et  al.  (2008b)  and  Parks  et  al.  (2014)  report  on  the 
2MASS  observations  of  a  small  region  in  Rho  Oph,  finding 
about  1 00  variables  with  roughly  similar  variability  properties  as 
Orion  in  that  the  amplitude  variations  were  found  to  be  between 
a  few  tenths  and  2  mag,  with  periods  obtained  for  about  one- 
third  of  the  sample.  Subsequent  to  2MASS,  more  recently,  there 
have  been  a  number  of  NIR  monitoring  programs  studying  other 
SFRs,  such  as  Wolk  et  al.  (2013),  which  monitored  Cyg  OB7, 
finding  several  classes  of  YSO  NIR  variability. 

In  the  mid-infrared  (Cohen  &  Schwartz  1 976),  as  for  the  near- 
infrared  (e.g.,  Elias  et  al.  1978;  Rydgren  &  Vrba  1983),  previous 
literature  suggested  at  least  small  variations  on  timescales 
of  months  to  years,  likely  attributable  to  circumstellar  disk 
processes.  However,  in  the  same  way  that  charge-coupled 
devices  (CCDs)  revolutionized  our  ability  to  discern  precisely 
optical  variability  trends  and  2MASS  did  the  same  for  near- 
infrared  variability,  the  Spitzer  Space  Telescope  (Werner  et  al. 
2004)  has  allowed  us  to  probe  even  small  variations  in  the  mid- 
infrared;  Spitzer  is  a  photometrically  stable  (better  than  1%), 
sensitive,  wide-field  (5'  x  5'),  Earth-trailing  (avoiding  orbital 
day /night  aliasing)  platform.  Spitzer,  specifically  the  Infrared 
Array  Camera  (IRAC;  Fazio  et  al.  2004),  observes  at  bands 
sensitive  to  both  YSO  photospheres  and  circumstellar  dust. 

Cycle  6  was  the  first  post-cryogen  Spitzer  cycle,  using  only 
IRAC’s  first  two  channels  (3.6  and  4.5  pm,  often  abbreviated 
as  IRAC-1  and  IRAC-2,  or  11  and  12;  when  reporting  mea¬ 
surements  in  magnitudes,  the  bands  are  written  with  brackets, 
e.g.,  [3.6]  =  16.38  mag).  The  YSOVAR  (Young  Stellar  Object 
VARiability)  the  Spitzer  Space  Telescope  Cycle-6  Exploration 
Science  (ES)  Program  was  approved  for  550  hr  of  observations, 
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Table  1 

Summary  of  Cluster  Properties 


Cluster a 

Dist. 

(pc) 

Gal.  Lat.b 

Class  II/I 
(G09)c 

Class  II/I 
(Obj.  w/L.C.)d 

Class  I/tot 
(Mem.  w/L.C.)e 

Notes 

AFGL  490 

900 

+  1.8 

-3.2 

4.5  ±  0.9 

0.20  ±  0.04 

Distance:  Testi  et  al.  (1998);  ratio:  Masiunas  et  al. 
(2012),  using  the  same  approach  as  G09  but  deeper 
data,  report  a  ratio  of  —5 

NGC  1333 

235 

-20:5 

-2.7 

4.7  ±  1.3 

0.16  ±0.04 

Distance:  Hirota  et  al.  (2008);  see  also  Hirota  et  al. 
(2011) 

Orion 

414 

—  19:0 

9.7  ±  0.9 

0.07  ±0.01 

Distance:  Menten  et  al.  (2007) 

Mon  R2 

830 

-12.6 

-4.7 

6.0  ±  1.4 

0.15  ±0.03 

Distance:  Herbst  &  Racine  (1976);  see  discussion  in 
Carpenter  &  Hodapp  (2008) 

GGD  12-15 

830 

-11.9 

-4.2 

5.8  ±  1.6 

0.13  ±0.03 

Distance:  Herbst  &  Racine  (1976);  see  discussion  in 
Carpenter  &  Hodapp  (2008) 

NGC  2264 

760 

+2.1 

3.0  ±  0.6 

0.15  ±0.03 

Distance:  Sung  et  al.  (1997) 

L1688 

120 

+16.6 

-3.0 

1.9  ±0.6 

0.28  ±  0.08 

Distance:  Wilking  et  al.  (2008),  Loinard  et  al. 

(2013),  Loinard  et  al.  (2008);  G09  cites  Wilking 
et  al.  (2005)  for  150  pc 

Serpens  Main 

415 

+16.5 

-1.4 

2.2  ±0.6 

0.21  ±0.05 

Distance:  Dzib  et  al.  (2010),  Loinard  et  al.  (2013) 
(see  also  Eiroa  et  al.  2008);  G09,  as  do  many  other 
authors,  cites  260  pc  from  Straizys  et  al.  (1996) 

Serpens  South 

415 

+3.8 

-0.7 

0.9  ±  0.2 

0.41  ±0.08 

G09  ratio  calculated  from  numbers  in  Gutermuth 
et  al.  (2008b);  assumed  to  be  same  distance  as 
Serpens  Main 

IRAS  20050+2720 

700 

-2.6 

-1.9 

2.2  ±  0.4 

0.30  ±  0.05 

Distance:  Wilking  et  al.  (1989),  Gunther  etal.  (2012) 

1C  1396A 

900 

+3.9 

7.9  ±  2.8 

0.04  ±0.01 

Distance:  Contreras  et  al.  (2002);  II/I  ratio  given  is 
in  region  with  both  11  and  12  Light  curves — II/I 
ratio  calculated  in  exactly  the  same  way  as  other 
clusters  is  1 1.6  ±  4.0 

Ceph  C 

700 

+2.1 

-2.3 

3.2  ±  1.0 

0.17  ±0.05 

Distance:  Moscadelli  et  al.  (2009);  II/I  ratio  given  is 
in  region  with  both  11  and  12  light  curves — II/I  ratio 
calculated  in  exactly  the  same  way  as  other  clusters 
is  3.8  ±  1.1 

Notes. 

a  Clusters  appear  in  R.A.  order,  in  this  and  subsequent  tables.  R.A.  for  our  specific  observations  appears  in  Table  2.  All  of  these  clusters  are  thought  to  be  between  1 
and  5  Myrold. 

b  Approximate  Galactic  latitude,  provided  as  a  rough  proxy  of  the  total  number  of  background/foreground  stars  (from  Galactic  contamination)  expected  in  the  region. 
c  Ratio  of  Class  II  to  Class  I  sources,  as  presented  in  G09.  For  a  brief  definition  of  SED  classes,  see  Appendix  B. 

d  Ratio  of  Class  11  to  Class  I  sources,  using  the  same  classification  scheme  as  G09  (where  the  classes  are  only  available  for  IR-selected  members),  but  using  the 
reprocessed  cryogenic  data,  and  the  ratio  is  calculated  only  for  objects  with  light  curves  in  the  YSOVAR  data.  Error  bars  are  derived  assuming  independent  Poisson 
uncertainties  on  the  numbers  used  to  compute  the  ratio.  See  Section  4.1  for  further  discussion. 

e  Ratio  of  Class  I  to  the  total  number  of  member  sources  with  viable  light  curves  in  the  YSOVAR  data,  using  the  standard  set  of  members  (Section  3.1)  and  our  SED 
classification  (Appendix  B).  Error  bars  are  derived  assuming  independent  Poisson  uncertainties  on  the  numbers  used  to  compute  the  ratio.  See  Section  4.1  for  further 
discussion  on  this  ratio  and  how  it  compares  to  the  Class  11/Class  I  ratio. 


with  the  goal  of  obtaining  the  first  extensive  mid-infrared  time 
series  photometry  of  the  central  ~1:  of  the  ONC  plus  smaller 
footprints  in  1 1  other  star-forming  cores;  see  Table  1  for  a  list 
of  the  clusters.  There  are  several  other  ES  and  smaller  programs 
exploring  YSOs  in  the  time  domain  with  Spitz.er,  many  of  which 
are  affiliated  (to  varying  degrees)  with  the  larger  YSOVAR  ef¬ 
fort.  We  have  incorporated  under  the  YSOVAR  umbrella  some 
additional  strongly  related  programs  pre-dating  and  arising  from 
YSOVAR,  resulting  in  a  total  of  786  hr  of  Spitzer  time.  About 
130  hr  of  that  is  dedicated  observations  of  NGC  2264  (Coordi¬ 
nated  Synoptic  Investigation:  NGC  2264,  or  CSI  2264),  which 
is  discussed  by  Cody  et  al.  (2014)  and  others  (e.g.,  Stauffer 
et  al.  2014;  Stauffer  et  al.  2015,  in  preparation).  CSI  2264  is 
not  discussed  in  the  same  way  as  the  rest  of  the  data  here,  in  no 
small  part  because  the  observations  are  generally  different  and 
because  it  involves  the  coordination  of  additional  telescopes.  A 
list  of  core  YSOVAR  programs  and  affiliated  programs  is  pre¬ 
sented  in  Table  2  and  discussed  in  Section  2.2.  We  sometimes 
refer  to  the  components  of  the  original  YSOVAR  program  as 
“YSOVAR-c lassie”  to  distinguish  them  from  the  smaller  affil¬ 


iated  observations  obtained  over  the  same  time  period.  There 
are  ~29,000  unique  objects  with  light  curves  from  either  (or 
both)  of  the  IRAC  channels  in  the  YSOVAR  data  set,  which  are 
matched  to  ~39,000  individual  light  curves.  These  light  curve 
counts  include  light  curves  from  both  cluster  members  and  a 
significant  number  of  non-member  stars,  and  also  likely  include 
extragalactic  objects.  There  are  more  light  curves  than  objects 
because  most  objects  have  light  curves  at  only  1 1  or  12,  but  many 
have  light  curves  at  both  II  and  12. 

YSOVAR  data  were  obtained  to  help  reveal  the  structure 
of  the  inner  disk  region  of  YSOs,  provide  new  constraints  on 
accretion  and  extinction  variability,  assess  timescales  of  mid-IR 
variability  from  seconds  to  years,  identify  new  young  eclipsing 
binaries,  help  identify  new  very  low-mass  substellar  members  of 
the  surveyed  clusters,  constrain  the  short-  and  long-term  stability 
of  hot  spots  on  the  surfaces  of  YSOs,  and  determine  rotational 
periods  for  objects  too  embedded  for  such  monitoring  in 
the  optical. 

In  this  paper,  in  addition  to  presenting  an  overview  of  the 
data  set,  one  of  our  goals  is  to  specifically  address  the  longest 
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Table  2 

Summary  of  YSOVAR  (and  Closely  Related)  Spitzer  Observations 


Cluster 

PIDa 

pib 

Approx.  Obs. 
Center  (J2000) 

Approx. 

Ecl.Lat. 

Date  Observed 

Epochs 

Cadeneec 

No.  Obj 
w/LCsd 

Total  Exp. 
Time  (hr) 

Notes 

AFGL  490 

60014 

J.  R.  Stauffer 

03:27:24  +58:44:00 

+38 

2011  Oct-Nov 

46 

Fast 

-1970 

9.1 

2  IRAC  FOVs 

NGC  1333 

61026 

J.  R.  Stauffer 

03:29:06  +31:19:30 

+  12 

2011  Oct-Nov 

73 

Fast 

-690 

18.5 

2x2  IRAC  FOVs 

Orion 

61028 

J.  R.  Stauffer 

05:35:15  -05:21:00 

-28 

2009  Oct-Dee; 

2010  Oct-Dee 

80 

Fast 

-7500 

365.2 

very  large  IRAC  map  (—0.9 
sq.  deg.) 

Orion 

61028 

J.  R.  Stauffer 

05:35:27  -04:47:31 

-28 

2011  Nov 

-270 

Staring 

-180 

10 

staring  eclipsing  binary  follow-up 

Orion 

61028 

J  R.  Stauffer 

05:35:27  -04:47:31 

-28 

2011  Nov 

94 

EB 

-480 

14.7 

2  FOVs,  mapping;  EB  follow-up 

Orion 

70025 

J.  R.  Stauffer 

05:35:02  -05:18:30 

-28 

2010  Nov-Dee 

74 

—Hourly 

e 

30.3 

Dipper  follow-up 

Mon  R2 

61025 

J.  R.  Stauffer 

06:07:48  -06:25:00 

-30 

2010  Nov-Dee 

46 

Fast 

-710 

5.6 

1  IRAC  FOV 

GGD  12-15 

61021 

J.  R.  Stauffer 

06:10:48  -06:12:30 

-30 

2010  Nov-Dee 

77 

Fast 

-1010 

14.8 

2  IRAC  FOVs 

GGD  12-15 

70172 

J.  Forbrich 

06:10:48  -06:12:30 

-30 

2010  Dec 

-530 

Staring 

-370 

20.0 

1  staring  FOV;  simul.  w/CXO,  & 
YSOVAR-likc  2-ficld  obs.  at 
start/end 

NGC  2264 

61027 

J.  R.  Stauffer 

06:41:04  +09:35:10 

-13 

2010  Nov-Dee 

39 

Fast 

-780 

4.7 

one  IRAC  FOV 

NGC  2264 

80040 

J.  R.  Stauffer 

06:40:48  +09:42:00 

-13 

2011  Dec 

f 

CSI 

-16,500 

99.1 

Cy8-CSI  2264,  staring  & 
mapping 

NGC  2264 

90098 

J  R.  Stauffer 

06:40:48  +09:42:00 

-13 

2013  Dec-2014 

Jan 

80 

CSI 

-4700 

30.3 

Cy9  -  CSI  2264  followup  on  -15 
targets  (cluster  targets) 

LI  688 

61024 

J  R.  Stauffer 

16:27:10-24:37:30 

-3 

2010  Apr-May; 

2010  Sep-Oct; 

201 1  Apr-May; 
2011  Oct-Nov 

108 

Fast/slow 

-840 

30.7 

3  IRAC  FOVs,  fast  cadence  2010 
Apr-May 

L1688 

60109 

P.  Plavchan 

16:27:31  -24:40:45 

-3 

2010  Apr 

-2500 

Staring 

-90 

24.0 

reverberation  mapping 

L1688 

90128 

H.  M.  Gunther 

16:27:31  -24:40:45 

-3 

2013  May-Jun 

10 

3-4  days 

-840 

2.8 

Cy9  follow-up 

Serpens  Main 

30319 

G.  Fazio 

18:29:59  +01:13:53 

+24 

2006  Sep 

29 

Cryo 

-2800 

6.5 

cryo  obs;  same  map  as  PC;  9  more 
hrs  staring  obtained  on  subset 

Serpens  Main 

61029 

J  R.  Stauffer 

18:29:59  +01:13:53 

+24 

201 1  May-Jun 

82 

Fasl 

-3400 

16.2 

2  IRAC  FOVs 

Serpens  South 

61030 

J.  R.  Stauffer 

18:30:04  -02:02:05 

+21 

201 1  May-Jun 

82 

Fast 

-1540 

10.0 

1  IRAC  FOV 

IRAS  20050+2720 

61023 

J.  R.  Stauffer 

20:07:04  +27:29:14 

+46 

2010  Jun-Aug 

102 

Fast 

-3030 

13.4 

1  IRAC  FOV 

IC  1396 A 

470 

B.  T.  Soifcr 

21:36:30  +57:29:48 

+64 

2008  Jan-Fcb 

29 

Cryo 

-4500 

15.0 

2x3  IRAC  FOVs;  DDT;  sec 
Morales-Calderon  ct  al.  (2009) 

IC  1396A 

497 

J.  R.  Stauffer 

21:36:45  +57:30:20 

+64 

2008  Dec 

10 

Cryo 

-1800 

1.5 

1  IRAC  FOV;  DDT;  16 
—contemporaneous  MIPS  24 pm 
epochs  obtained,  2.1  more  hrs. 

IC  1396 A 

61022 

J.  R.  Stauffer 

21:36:30  +57:29:48 

+64 

2009  Aug-2010 
Mar;  2010 
Aug-2011  Feb 

143s 

Fast/slow 

-5100 

38.9 

2x3  IRAC  FOVs;  fast  cadence 
2009  Scp-Nov 

Ccph  C 

61020 

J.  R.  Stauffer 

23:05:51  +62:30:55 

+59 

2009  Dec-2010 
Mar;  2011 
Jan-Mar;  201 1 
Scp-Nov 

39 

Fast/slow 

-1950 

4.8 

1  IRAC  FOV;  fast  cadence  2010 
Aug-Nov 

(Ccph  C) 

(61020) 

(K.  Covey) 

23:05:51  +62:30:55 

+59 

2010  Aug-Nov 

105 

Fast 

(as  above) 

13.5 

1  FOV;  simul.  w/CXO  (CXO 

Cyl  1  PI  :K.  Covey);  AORs 
entirely  included  within  program 
61020 

Notes. 

a  PID  =  Program  IDcntification  number;  data  from  each  set  of  observations  of  each  cluster  was  grouped  within  the  same  program.  The  data  associated  with  each  PID  can  be  retrieved  all  at 
once  from  the  Spitzer  Heritage  Archive  using  this  number.  YSOVAR-classic  programs  arc  60014  and  61020-61030,  inclusive. 
b  PI  =  Program’s  Primary  Investigator. 

c  Cadence,  meaning  rate  at  which  observations  were  obtained.  For  the  YSOVAR-classic  clusters,  the  fast/slow  cadence  is  described  in  Section  2.4.  Additional  values  found  in  this  column 
include  staring  (meaning  that  the  IRAC  observations  were  staring  rather  than  mapping  mode),  EB  (observations  designed  for  follow  up  of  specific  eclipsing  binary  candidates),  hourly,  CS1 
(cadence  described  in  Cody  ct  al.  2014),  cryo  (observations  conducted  in  cryogenic  era  and  cadence  was  different  for  each  program). 

d  Approximate  number  of  unique  objects  with  a  light  curve  in  cither  or  both  IRAC  channels.  Counting  only  the  original  YSOVAR  programs,  there  arc  —29,000  unique  objects  with  a  light 
curve  in  cither  or  both  IRAC  channels.  Counting  IRAC- 1  and  IRAC-2  for  the  same  object  as  two  distinct  light  curves,  there  arc  —  1 1 ,000  objects  with  light  curves  in  both  channels,  and  a  total 
of  —39,000  light  curves. 

c  Counted  as  part  of  light  curves  listed  for  program  61028. 
f  344  AORs  used  in  this  program. 

s  One  more  epoch  was  also  obtained  during  the  IRAC  Warm  Instrument  Calibration  (IWIC),  which  was  before  IRAC  was  properly  calibrated  during  the  initial  days  of  the  warm  mission,  and 
so  that  epoch  is  not  included  in  our  final  data  set. 


timescale  variations  that  we  can  quantify  in  these  clusters,  six 
to  seven  years.  We  look  for  large  changes  between  the  earliest 
Spitzer  observations  (from  observations  obtained  early  in  the 
cryogenic  era)  and  the  YSOVAR  monitoring  observations.  We 
discuss  many  of  the  technical  details  associated  with  the  data 
reduction  across  the  YSOVAR  effort.  The  Orion  data  were  first 
described  by  Morales-Calderon  et  al.  (201 1,  hereafter  MCI  1). 
The  other  1 1  smaller-field  clusters  are  introduced  in  this  paper, 
but  will  be  discussed  in  detail  in  other  papers  (the  first  of  which, 
on  LI 688,  is  Gunther  et  al.  2014).  Here,  we  first  provide  a 


summary  of  the  observations  and  data  reduction  (Section  2), 
followed  by  a  definition  of  the  samples  we  use  (Section  3).  We 
present  some  global  statistics  on  all  the  clusters  in  Section  4.  We 
delve  into  variable  selection  in  Section  5,  discussing  different 
tests  for  selecting  variables  and  simulating  the  sensitivity  of 
the  techniques  given  the  actual  data  set.  We  present  some 
analysis  that  is  best  performed  with  all  the  clusters  together 
in  Section  6,  such  as  fractions  of  long-term  variables  across 
clusters,  correlations  between  rotation  rate  and  IR  excess,  and 
the  absence  of  transient  disks.  We  summarize  in  Section  7. 
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For  completeness,  we  note  here  that  we  refer  to  the  12  re¬ 
gions  of  recent  star  formation  that  we  observed  for  YSOVAR 
as  “clusters,”  knowing  that  others  may  prefer  “associations” 
or  other  nomenclature.  Our  targets  resemble  small  condensa¬ 
tions  within  a  region;  Gutermuth  et  al.  (2009),  among  others, 
addresses  formal  clustering  in  these  regions. 

2.  OBSERVATIONS  AND  DATA  REDUCTION 

In  this  section,  we  review  the  target  selection  and  some 
general  properties  of  the  clusters,  and  describe  the  Warm 
Spilzer  observations:  observing  strategy,  data  reduction,  ca¬ 
dence,  and  the  noise  floor.  We  also  describe  the  data  reduction 
for  the  cryogenic-era  data,  other  archival  data,  and  the  Chandra 
X-ray  data. 

2.1.  Target  Selection 

In  this  section,  we  first  discuss  our  criteria  for  picking  targets 
and  then  review  basic  properties  of  each  cluster. 

2.1.1.  Overview  of  Target  Selection 

The  Orion  SFR  has  been  the  subject  of  variability  studies 
more  than  any  other  SFR.  Because  Orion  is  so  well-studied, 
particularly  for  optical  and  NIR  variability,  we  elected  to  include 
it  as  a  very  significant  part  of  our  mid-infrared  observing  effort. 

Besides  Orion,  we  selected  the  cores  of  1 1  additional  young 
clusters  for  monitoring;  see  Table  1.  Our  smaller-field,  more 
embedded  cluster  sample  (“smaller  field  clusters”)  was  selected 
based  on  detailed  examination  of  all  the  SFRs  surveyed  with 
Spitz.er  by  mid-2007,  when  YSOVAR  targets  were  selected. 
We  chose  regions  that  satisfy  the  following  criteria.  (1)  A 
relatively  high  fraction  of  Class  I  sources  (for  a  brief  definition 
of  spectral  energy  distribution  (SED)  classes,  such  as  Class  I, 
see  Appendix  B),  such  that  we  would  obtain  monitoring  for 
some  of  the  most  heavily  embedded  objects;  (2)  a  high  density 
of  YSOs  within  one  or  a  few  IRAC  fields  of  view  (FOVs)  - 
typically  >40  YSOs  per  field;  (3)  moderate  cirrus  backgrounds; 
and  (4)  minimal  problems  with  crowding  or  very  bright  nearby 
sources.  Several  of  the  regions  we  monitored  are  very  difficult  to 
observe  from  the  ground  due  to  their  high  level  of  obscuration. 
We  included  IC  1396A  even  though  it  is  not  as  embedded 
because  there  has  already  been  a  significant  investment  of 
Spitz, er  time  in  monitoring  this  region  (more  details  on  this 
follow  below).  NGC  2264,  like  Orion,  has  been  intensively 
surveyed  for  variability  from  the  ground  (see,  e.g.,  Makidon 
et  al.  2004;  Lamm  et  al.  2003;  Cieza  &  Baliber  2007)  and  from 
space  (see,  e.g.,  Alencar  et  al.  2010  for  CoRoT;  Zwintz  et  al. 
2009  for  MOST).  We  focused  on  the  most  embedded  region  of 
NGC  2264  for  the  “YSOVAR-classic”  part  of  the  program  and 
monitored  a  much  larger  region  for  CSI  2264. 

2.1.2.  Cluster  Properties 

We  now  briefly  discuss  general  properties  of  each  of  these 
clusters,  in  R.A.  order  (see  Tables  1  and  2).  We  tabulate  several 
characteristics  of  these  clusters  in  Table  1 ,  including  some  values 
from  Gutermuth  et  al.  (2009,  2010,  hereafter  G09).  All  of  these 
clusters  are  thought  to  be  between  1  and  5  Myr  old.  Ages  more 
accurate  than  that  for  clusters  as  young  and  embedded  as  these 
are  difficult  to  obtain,  even  in  a  relative  sense.  Cryogenic  Spitzer 
observations  for  many  of  these  clusters  were  discussed  in  G09, 
who  calculated  the  Class  II  to  Class  I  ratio  for  various  regions 
and  sub-clusters  within  the  cryo  Spitzer  maps.  This  ratio  is  easier 


to  obtain  than  an  age  and  helps  place  the  clusters’  evolutionary 
states  into  context  with  each  other;  it  is  provided  here  (in  Table  1 
and  the  text  below)  in  part  as  a  link  back  to  existing  literature. 
In  general,  our  observations  enclose  only  the  most  embedded 
parts  of  these  clusters;  the  full  cluster  membership  generally 
includes  objects  beyond  the  regions  we  monitored.  (Orion  is 
different  in  that  the  map  is  much  larger  than  the  other  smaller- 
field  cluster  maps,  but,  even  then,  the  map  does  not  include  all 
objects  thought  to  be  part  of  Orion.)  Thus,  we  have  recalculated 
the  Class  II  to  Class  I  ratio  for  only  objects  with  YSOVAR  light 
curves  in  the  same  way  as  was  done  in  G09  (e.g.,  only  for  the 
IR-selected  members);  this  value  appears  in  Table  1,  again  to 
place  our  observations  into  context  with  the  prior  literature.  We 
discuss  this  parameterization  in  more  detail  below  in  Section  4. 1 , 
and  consider  a  different  parameterization  of  the  clusters’  relative 
evolutionary  states,  the  ratio  of  the  number  of  Class  I  sources 
to  total  YSOs;  anticipating  this  discussion,  these  values  are 
included  in  Table  1 .  The  total  cluster  membership  we  use  for 
this  metric  includes  objects  with  light  curves  that  we  selected 
using  both  the  IR  data  from  Spitzer  and  X-ray  data  from  the 
Chandra  X-ray  Observatory  (Section  3).  The  X-ray  data  are 
discussed  further  in  Section  2. 10. 

Table  1  also  includes  an  approximate  Galactic  latitude  as  a 
very  rough  guide  to  the  overall  Galactic  background/foreground 
source  density  expected  near  each  location.  Targets  close  to  the 
Galactic  plane  (e.g.,  Serpens  Main)  have  a  higher  surface  density 
of  objects  in  our  FOVs  than  targets  further  from  the  plane. 

AFGL  490.  The  cluster  associated  with  AFGL  490  is  a  portion 
of  the  larger  Cam  OB  1  Association,  centered  on  a  massive  object 
(8-10  MQ\  sometimes  the  name  AFGL  490  is  used  to  refer  only 
to  this  massive  protostar).  The  cluster  is  thought  to  lie  between 
~900  pc  (Testi  et  al.  1998)  and  ~  1 0 1 0  pc  (Straizys  &  Laugalys 
2008).  We  adopt  a  distance  of  900  pc  (as  do  G09).  G09,  based 
on  the  cryogenic  Spitzer  data  for  this  region,  found  at  least  100 
young  stars  or  candidates  in  this  vicinity,  and  found  the  overall 
Class  II  to  Class  I  ratio  to  be  about  3.  Masiunas  et  al.  (2012)  also 
report  on  the  cryogenic  Spitzer  observations  and,  incorporating 
additional  data,  find  many  new  YSO  candidates  and  a  Class  II 
to  Class  I  ratio  overall  of  ~5.  Among  the  objects  with  YSOVAR 
light  curves,  using  the  G09  color  cuts  to  identify  Class  I  and 
II  sources  (see  Section  2.7  and  Appendix  B),  the  ratio  is  ~4.5. 
This  region  is  the  only  cluster  in  our  set  of  clusters  that  does  not 
have  archival  Chandra  X-ray  observations  (see  Section  2.10). 

NGC  1333.  NGC  1333  is  on  the  western  edge  of  the  Perseus 
molecular  cloud,  and,  at  only  ~235  pc  (Hirota  et  al.  2008, 2011), 
it  is  one  of  the  youngest  and  most  well-studied  SFRs,  with 
>200  refereed  publications.  The  region  is  riddled  with  outflows 
from  young  stars  (e.g.,  Plunkett  et  al.  2013),  which  can  be  seen 
clearly  in  the  Spitzer  4.5  /am  image  of  this  region.  Gutermuth 
et  al.  (2008b,  2009)  analyzed  the  cryogenic  Spitzer  data  for  this 
region,  and  found  more  than  130  young  stars  or  candidates; 
the  Class  II  to  Class  I  ratio  they  found  in  the  core  is  ~2.7. 
For  the  region  where  we  have  light  curves,  the  ratio  is  ~4.7. 
Despite  the  large  number  of  previous  studies  in  the  literature, 
NGC  1333  has  few  prior  studies  specifically  investigating  time 
variability.  The  YSOVAR  data  will  be  discussed  in  detail  in 
L.  Rebull  et  al.  (in  preparation);  Raga  et  al.  (2013)  report 
on  proper  motions  of  the  outflows  in  this  region  using  the 
YSOVAR  data. 

Orion.  Orion  has  been  the  subject  of  extensive  variability 
studies.  Haro  (1969),  Herbig  &  Kameswara  Rao  (1972),  and 
Walker  (1978)  all  conducted  pioneering  studies  of  the  variability 
of  young  stars  in  Orion,  leading  to  the  modern  era  using 
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two-dimensional  imaging  cameras  as  initiated  by  Herbst  and 
his  team  (Attridge  &  Herbst  1992;  Choi  &  Herbst  1996). 
They  obtained  multi-year,  optical  time  series  photometry  of 
several  regions  of  the  ONC,  eventually  obtaining  light  curves 
for  hundreds  of  YSOs,  often  spanning  several  years.  Many 
other  groups  subsequently  obtained  time  series  photometry  of 
other  portions  of  the  ONC,  using  optical  (Stassun  et  al.  1999, 
2006,  2007;  Rebull  2001;  Herbst  et  al.  2002;  Irwin  et  al.  2007; 
Rodrfguez-Ledesma  et  al.  2009)  and  near-IR  (CHS01)  imaging 
data.  More  recently,  some  far-IR  monitoring  has  been  conducted 
in  Orion  as  well  (Billot  et  al.  2012).  In  many  cases,  the  light 
curve  shapes  are  well-fitted  by  models  of  rotational  modulation 
via  hot  or  cold  spots,  normally  at  moderately  high  latitudes 
since  the  light  curves  are  seldom  “flat-bottomed.”  However,  for 
a  significant  fraction  of  the  light  curves,  particularly  those  in  the 
near-IR,  spots  do  not  seem  to  provide  a  good  explanation  for 
the  observed  variability  (CHS01).  We  included  the  ONC  as  a 
YSOVAR  target  because  of  the  substantial  amount  of  extant 
monitoring  available  in  the  literature.  Note  that  it  is  likely 
slightly  older  than  most  of  the  other  embedded  regions  studied 
here.  For  the  cryogenic-era  data  (see  Section  2.7),  we  used  the 
data  reduction,  YSO  identification,  and  YSO  classification  from 
Megeath  et  al.  (2012),  which  are  very  similar  to  that  from  G09. 
Therefore,  while  there  is  no  ratio  of  Class  II  to  Class  I  objects 
from  G09  to  report,  we  calculated  this  ratio  using  the  Megeath 
et  al.  (2012)  classifications  for  only  those  objects  for  which  we 
have  light  curves,  obtaining  ~9.7.  While  this  ratio  is  affected 
by  the  much  larger  region  monitored  by  YSOVAR  (only  the 
cores  are  monitored  in  most  of  the  other  regions),  this  value  is 
consistent  with  Orion  being  slightly  older  than  our  other  more 
embedded  targets.  We  adopt  a  distance  of  414  pc  from  Menten 
et  al.  (2007). 

Mon  R2.  Mon  R2  appears  in  G09,  with  a  Class  II  to  Class  I 
ratio  in  the  central  region  of  ~4.7.  This  region,  part  of  the 
Monoceros  R2  molecular  cloud,  is  near  vdB  67  and  69  (see, 
e.g..  Carpenter  &  Hodapp  2008).  It  is  typically  thought  to  be  at 
~830  pc  (Herbst  &  Racine  1976;  Carpenter  &  Hodapp  2008). 
There  are  several  sources  that  are  very  bright  in  the  infrared 
in  this  location,  which  affect  the  completeness  of  the  catalog 
extracted  from  the  Spitzer  data.  The  Class  II  to  Class  I  ratio  in 
the  region  with  YSOVAR  monitoring  is  ~6.  These  data  will 
be  discussed  in  more  detail  by  L.  A.  Hillenbrand  et  al.  (in 
preparation). 

GGD  12-15.  GGD  12-15  is  a  dense  core  also  located  in  the 
Monoceros  R2  molecular  cloud  (see,  e.g.,  Carpenter  &  Hodapp 
2008).  As  such,  we  assume  it  to  also  be  at  ~830  pc.  G09  find  a 
Class  II  to  Class  I  ratio  of  ~4.2;  we  calculate  ~5.8  for  the  region 
with  YSOVAR  light  curves.  Several  time-variable  radio  sources 
are  located  here  (Carpenter  &  Hodapp  2008,  and  references 
therein).  We  also  monitored  this  region  with  Chandra  during 
portions  of  our  YSOVAR  Spitzer  campaign.  These  data  will  be 
discussed  in  depth  by  S.  Wolk  et  al.  (in  preparation). 

NGC  2264.  NGC  2264  is  thought  to  be  comparable  in  age  to 
or  slightly  older  than  Orion  (see,  e.g.,  Ramirez  et  al.  2004a).  This 
region  does  not  appear  in  G09  so  we  do  not  have  a  similarly 
obtained  Class  II  to  Class  I  ratio  for  comparison.  However, 
the  region  we  monitored  as  part  of  the  YSOVAR-classic  data 
(original  YSOVAR  program  6 1 027 ;  see  Section  2.2  and  Table  2), 
which  is  the  region  discussed  here,  is  centered  on  the  Spokes 
Cluster  (Teixeira  et  al.  2006),  the  most  embedded  portion  of 
NGC  2264  (analogous  to  the  BN  region  in  Orion).  We  calculate 
a  Class  II  to  Class  I  ratio  in  this  monitored  region  of  ~3.0, 
which  is  comparable  to  the  ratio  for  many  of  the  other  embedded 


clusters  here.  We  work  only  with  the  YSOVAR-classic  data  in 
the  present  paper;  Cody  et  al.  (2014)  discuss  the  much  larger 
CSI 2264  data  set.  We  have  adopted  a  distance  to  this  cluster  of 
760  pc  (Park  et  al.  2000). 

LI 688.  Lynds  1688  (LI 688)  is  located  within  the  p  Ophiuchi 
molecular  cloud.  It  is  also  one  of  the  best-studied  SFRs  in  this 
YSOVAR  data  set,  with  more  than  500  refereed  articles.  The 
distance  to  this  region  is  a  subject  of  some  debate.  There  is 
recent  evidence  for  130  pc  (Wilking  et  al.  2008,  and  references 
therein),  131  pc  (Mamajek  2008),  and  120  pc  (Loinard  et  al. 
2008,  2013;  Lombardi  et  al.  2008).  Lombardi  et  al.  (2008) 
offer  a  plausible  explanation  for  the  “discrepancy”  in  the  very 
long  baseline  interferometry  (VLBI)  results  noted  by  Wilking 
et  al.  (2008);  there  may  be  other  subregions  in  Oph  at  distinct 
distances,  but  the  evidence  is  unclear  at  this  point.  Our  results  are 
not  particularly  dependent  on  distance;  we  have  adopted  120  pc 
as  one  of  the  most  recent  determinations.  Lynds  1688  has  a 
high  surface  density  of  embedded  objects.  This  region  appears 
in  G09  with  a  Class  II  to  Class  I  ratio  of  ~3.0;  we  calculate 
a  value  of  ~1.9  for  the  region  we  monitored.  Barsony  et  al. 
(2005)  first  reported  YSO  variability  in  the  MIR  in  the  objects 
in  this  core.  By  comparison  to  Infrared  Space  Observatory 
data,  they  found  significant  variability  in  18  out  of  85  objects 
detected,  on  timescales  of  years.  They  found  such  variability  in 
all  SED  classes  with  optically  thick  disks,  and  suggest  that  this 
might  be  due  to  time-variable  accretion.  The  Multiband  Imaging 
Photometer  for  Spitzer  (MIPS;  Rieke  et  al.  2004)  data  for  this 
region  were  presented  in  Padgett  et  al.  (2008);  no  variability  in 
sources  at  24  pm  was  found  to  a  level  of  10%  on  timescales 
of  hours.  Alves  de  Oliveira  &  Casali  (2008)  recently  reported 
on  deep  NIR  monitoring  of  this  region,  with  14  epochs  over 
two  years,  finding  that  41%  of  the  known  YSOs  are  variable. 
The  p  Oph  core  region  was  included  in  one  of  the  “calibration” 
2MASS  fields,  and  as  such  has  monitoring  data  in  the  NIR  (Parks 
et  al.  2014)  with  a  cadence  of  ~1  day  over  3  observing  seasons 
spanning  ~2.5  yr.  Parks  et  al.  (2014)  found  101  variables,  72  of 
which  are  identified  with  known  YSOs  in  the  region.  Plavchan 
et  al.  (2013)  report  on  YLW  16A,  finding  a  93  day  periodicity. 
The  YSOVAR  data  are  discussed  in  detail  by  Gunther  et  al. 
(2014). 

Serpens  Main.  The  Serpens  core  has  been  studied  for  decades, 
but  has  become  known  as  “Serpens  Main”  to  distinguish  it 
from  the  relatively  recently  discovered  embedded  star-forming 
core  known  as  Serpens  South  (Gutermuth  et  al.  2008b;  see 
below).  The  original  Spitzer  cry  o-era  data  for  Serpens  Main  were 
presented  in  Harvey  et  al.  (2007);  they  found  no  clear  evidence 
at  the  ~25%  level  for  IRAC-band  variability  in  any  sources  in 
the  field  over  the  ~6  hr  timescale  of  their  observations.  The 
Spitzer  cryogenic  data  were  also  used  in  G09,  who  determined 
the  ratio  of  Class  II  to  Class  I  objects  at  ~1.4,  the  lowest  of  all 
of  the  clusters  from  YSOVAR  that  also  appear  in  G09.  It  has 
a  high  surface  density  of  embedded  objects  and  it  is  another 
very  well-studied  SFR,  with  more  than  500  refereed  articles. 
While  the  Straizys  et  al.  (1996)  distance  of  260  pc  to  Serpens 
Main  is  well-cited  in  the  literature,  more  recent  studies  (Dzib 
et  al.  2010;  Loinard  et  al.  2013;  see  also  Eiroa  et  al.  2008) 
using,  e.g.,  VLBI  instead  suggest  that  the  distance  of  260  pc 
may  be  portions  of  clouds  associated  with  Aquila,  and  that 
a  better  distance  for  Serpens  itself  is  actually  415  pc,  which 
we  adopt  here.  Hodapp  (1999)  reports  on  NIR  variability  of 
knots,  jets,  and  young  stars;  in  terms  of  point  sources,  Hodapp 
(1999)  primarily  discusses  one  particular  source  (OO  Ser)  in 
this  region.  There  was  a  brief  Spitzer  monitoring  program  of 
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this  region  conducted  during  Cycle  3  as  part  of  guaranteed  time 
(see  Table  2);  the  observations  conducted  as  part  of  YSOVAR 
were  designed  to  be  well-matched  to  these  observations.  We 
obtain  a  Class  II  to  Class  I  ratio  of  —2.2  for  the  region  with  light 
curves. 

Serpens  South.  Serpens  South  was  discovered  by  Gutermuth 
et  al.  (2008b)  as  a  dense,  embedded  cluster  in  the  Serpens- 
Aquila  Rift.  It  is  thought  that  this  cluster  is  at  about  the  same 
distance  as  Serpens  Main,  which  we  have  taken  to  be  415  pc. 
From  the  numbers  in  Gutermuth  et  al.  (2008b),  it  has  a  Class  II 
to  Class  I  ratio  of  only  ~0.7.  Considering  only  the  region  we 
monitored,  the  ratio  is  comparable  at  ~0.9. 

IRAS  20050+2720.  Observations  from  Spitzer  and  Chandra 
of  the  cluster  associated  with  IRAS  20050+2720  (abbreviated 
IRAS  20050)  have  been  discussed  by  Gunther  et  al.  (2012). 
This  cluster  is  part  of  the  Cygnus  Rift  and  is  likely  at  ~700  pc 
(Gunther  et  al.  2012,  and  references  therein).  G09  determined 
the  ratio  of  Class  II  to  Class  I  objects  to  be  — 1.9;  in  the  region 
with  light  curves,  we  obtain  ~2.2.  The  YSOVAR  data  for  this 
cluster  will  be  discussed  in  detail  by  K.  Poppenhaeger  et  al.  (in 
preparation). 

IC  1396A.  IC  1396 A  is  the  most  prominent  globule  in  the 
IC  1396  complex;  this  cluster  is  sometimes  called  the  Elephant 
Trunk  Nebula.  Morales-Calderon  et  al.  (2009)  presented  the 
first  Spitzer  monitoring  of  young  stars,  centered  on  this  tar¬ 
get.  This  region  is  also  likely  to  be  very  young,  but  it  is  not 
particularly  embedded,  at  least  compared  to  other  YSOVAR 
clusters.  It  is  not  included  in  G09,  but  over  the  entire  re¬ 
gion  within  which  light  curves  were  obtained,  we  obtain  a 
Class  II  to  Class  I  ratio  of  ~  1 1 .6.  Because  IC  1396A  is  at 
a  high  ecliptic  latitude  (see  Section  2.3  below  for  discus¬ 
sion  of  ecliptic  latitude  dependencies),  there  are  light  curves 
with  only  a  few  single-band  points  for  many  objects.  For  the 
much  smaller  number  of  objects  that  have  light  curves  in  both 
IRAC-1  and  -2,  we  obtain  a  much  lower  Class  II  to  Class  I  ratio 
of  ~7.9,  though  within  Poisson  errors  calculated  assuming  inde¬ 
pendent  errors  in  the  numerator  and  denominator,  these  values 
are  consistent  (1 1.6  ±  4.0  and  7.9  ±  2.8).  We  assume  it  is  at  a 
distance  of  —900  pc  (Contreras  et  al.  2002). 

Ceph  C.  Ceph  C  is  part  of  the  Cep  OB  3  molecular  cloud, 
and  is  included  in  the  G09  study,  with  a  Class  II  to  Class  I 
ratio  of  —2.3  obtained  in  that  paper.  G09  used  a  distance  of 
730  pc  (Blauw  1964);  we  adopt  a  distance  of  700  pc  based  on 
maser  parallax  from  Moscadelli  et  al.  (2009).  Ceph  C  is  one 
of  the  less  well-studied  clusters  in  our  set.  For  the  region  with 
any  light  curves  (see  Section  2.3  below),  we  obtain  a  Class  II 
to  Class  I  ratio  of  ~3.8.  Because  it,  like  IC  1396A,  is  at  a 
high  ecliptic  latitude,  we  repeated  this  calculation  for  the  much 
smaller  number  of  objects  that  have  light  curves  in  both  IRAC- 
1  and  -2,  obtaining  a  ratio  of  —3.2;  again  assuming  Poisson 
counting  statistics,  these  ratios  are  comparable  (3.8  ±  1.1  and 
3.2  ±  1.0).  This  cluster  was  monitored  during  our  YSOVAR 
campaign  with  Chandra  as  well;  these  data  will  be  discussed  in 
depth  by  K.  Covey  et  al.  (in  preparation). 

2.2.  Warm  Spitzer  Observations:  General  Properties 

To  better  manage  the  observation  planning  and  data  down¬ 
loading  for  the  12  clusters  we  observed,  we  separated  each 
cluster  into  an  individual  observing  program.  The  individual 
observing  programs  and  clusters  are  listed  in  Table  2;  one  can 
download  the  data  from  the  Spitzer  Heritage  Archive  using  these 
program  numbers.  We  will  deliver  all  of  our  extracted  photom¬ 
etry  to  the  Spitzer  Science  Center  (SSC)  and  Infrared  Science 


Archive  (IRS  A29)  for  public  access  via  IRS  A  tools,  and  through 
those  tools,  the  Virtual  Observatory. 

Most  of  the  observations  were  IRAC  mapping  mode  Astro¬ 
nomical  Observation  Requests  (AORs)  using  full-array  12  s 
high-dynamic-range  (HDR)  mode  (which  is  defined  such  that  a 
0.4  and  10.4  s  exposure  is  obtained  at  each  pointing — see  the 
IRAC  Instrument  Handbook,  available  at  the  SSC/IRSA  Web 
site30).  In  some  cases,  as  noted  in  Table  2,  the  AORs  were  staring 
AORs,  meaning  that  they  were  continuous  or  semi-continuous 
observations  without  dithering  or  mapping.  These  staring  data 
will  be  discussed  in  other  YSOVAR  papers.  The  present  paper  is 
limited  to  the  HDR  mapping  observations  (and  is  largely  further 
restricted  to  the  “fast  cadence”  observations;  see  Section  2.4). 

Roughly  half  of  the  original  YSOVAR  time  allocation  was 
devoted  to  observations  of  Orion.  The  Orion  Spitzer  campaign 
included  a  —0.9  deg2  region  centered  on  the  Trapezium  cluster 
in  Orion;  this  is  far  larger  than  can  be  obtained  in  a  single 
AOR  at  a  single  epoch.  The  observed  area  was  thus  broken 
into  five  segments  with  a  central  region  of  —20'  x  25'  and 
four  flanking  fields.  The  central  part  was  observed  in  full  array 
mode  with  1.2  s  of  exposure  time  and  20  dither  positions  to 
avoid  saturation  by  the  bright  nebulosity  around  the  Trapezium 
stars.  The  remaining  four  segments  of  the  map  were  observed 
in  12  s  HDR  mode,  and  four  dither  positions.  More  details  on 
these  original  YSOVAR  Orion  observations  appear  in  MC11. 
We  reallocated  time  within  the  YSOVAR  time  budget  to  follow 
up  on  some  eclipsing  binaries  in  Orion  (Morales-Calderon  et  al. 
2012).  We  also  obtained  time  in  Cycle  7  to  follow  up  on  some 
AA  Tau  analogues  in  Orion  presented  in  MC 1 1 .  These  follow-up 
observations  did  not  map  the  entire  Orion  region. 

The  rest  of  the  original  YSOVAR  time  allocation  was  divided 
among  the  rest  of  the  SFRs.  In  contrast  to  Orion,  for  the  1 1  other 
clusters,  the  monitoring  observations  are  one  or  at  most  a  few 
IRAC  FOVs;  see  the  last  column  of  Table  2.  As  noted  above, 
we  sometimes  refer  to  these  1 1  other  clusters  as  “the  smaller 
field  clusters.” 

The  observations  are  typically  spread  over  weeks  to  months, 
usually  at  a  cadence  of  about  twice  per  day,  but  with  an  interval 
between  observations  that  varied  both  by  design  (to  reduce 
aliasing  problems)  and  due  to  constraints  imposed  by  other 
Spitzer  programs  or  infrastructure  operations  that  were  executed 
during  the  same  campaigns  (further  details  on  the  cadence 
follow  below  in  Section  2.4). 

NGC  2264  was  included  as  a  smaller-field  cluster  to  be 
monitored  as  part  of  the  original  YSOVAR  program,  and 
while  this  original  program  was  still  executing,  CSI  2264  was 
approved.  We  continued  to  execute  the  original  YSOVAR  small- 
field  observations  in  this  region  (program  61027  in  Table  2),  and 
these  small-field  regions  are  discussed  here,  since  they  resemble 
the  rest  of  the  original  YSOVAR  programs  more  closely  than 
they  resemble  CSI  2264. 

Two  of  our  clusters,  IC  1396 A  and  Serpens  Main,  were 
monitored  in  the  cryogenic  era  with  Spitzer  with  the  primary 
intention  of  monitoring  changes  in  these  objects  in  the  Spitzer 
bands  (as  opposed  to  removing  artifacts).  In  YSOVAR,  we  re¬ 
observed  these  clusters  in  the  same  way  so  that  the  same  objects 
continue  to  be  monitored,  and  we  re-reduced  all  of  these  data 
in  the  same  fashion  as  discussed  here.  In  the  context  of  this 
paper,  we  are  not  particularly  focused  on  these  cryogenic-era 


29  http://irsa. ipac.caltech.edu/ 

30  http://irsa. ipac.caltech.edu/data/SPITZER/docs/ 
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Figure  1.  Approximate  sky  coverage  for  a  summed-up  image  consisting  of 
all  epochs  of  YSOVAR  AFGL  490  observations,  superimposed  on  a  reverse 
grayscale  image  of  AFGL  490  at  4.5  pm  obtained  during  the  cryogenic  mission. 
The  thicker  solid  blue  line  is  3.6  pm  and  the  thicker  red  dashed  line  is  4.5  pm. 
A  single  epoch  of  observation  is  also  indicated  by  thinner  solid  blue  and  dashed 
red  lines,  with  the  difference  between  the  single  epoch  and  the  larger  polygon 
due  to  (ecliptic-latitude-dependent)  field  rotation  effects.  North  is  up  and  east  is 
to  the  left.  The  distance  between  the  farthest  north  and  farthest  south  coverage 
here  is  ~27'.  The  field  rotation  here  is  not  as  significant  as  in  those  with  ecliptic 
latitudes  ~60  .  There  is  no  Chandra  coverage  for  this  cluster.  Similar  figures 
for  the  remaining  1 1  clusters  are  included  in  Appendix  C. 

(A  color  version  of  this  figure  is  available  in  the  online  journal.) 


monitoring  data;  however,  these  data  will  be  included  in  the 
YSOVAR  cluster-specific  papers. 

2.3.  Warm  Spitzer  Observations:  Footprints  and 
Operational  Constraints 

Because  of  the  nature  of  Spitzer  and  IRAC,  the  footprints 
of  our  observations  and  how  they  vary  with  time  are  both 
complicated  issues.  We  now  discuss  these  issues  as  they  pertain 
to  YSOVAR  observations. 

Figures  1-13  present  the  outline  (footprint)  of  the  YSOVAR 
observations.  The  original  YSOVAR  Orion  observations  (Fig¬ 
ure  3)  and  the  CSI 2264  observations  (Figure  7)  cover  a  substan¬ 
tially  larger  area  than  those  for  the  other  smaller-field  clusters. 
Chandra  footprints  are  included  in  these  figures  for  reference 
and  discussed  below  in  Section  2.10. 

The  focal  plane  of  the  IRAC  camera  is  such  that  the  3.6  and 
4.5  pm  FOVs  are  not  the  same;  data  are  obtained  in  both  FOVs 
at  once,  with  one  field  placed  on  the  target  of  interest,  and  the 


Figure  2.  Same  as  Figure  1,  but  for  our  NGC  1333  observations,  with  the 
north-south  boundary  extremes  separated  by  ~24'.  The  field  rotation  here  is 
essentially  zero.  The  yellow  polygon  indicates  the  approximate  region  covered 
by  Chandra  observations. 

(A  color  version  of  this  figure  is  available  in  the  online  journal.) 


other  obtaining  serendipitous  data  in  a  non-overlapping  ~5'  x  5' 
field  with  a  center  offset  ~6.'5  from  the  target  field;  see  the 
IRAC  Instrument  Handbook  for  more  details.  The  placement  of 
these  FOVs  also  changes  with  time.  For  all  Spitzer  observations, 
as  discussed  in  the  Spitzer  Space  Telescope  Handbook  (also 
available  at  the  SSC/IRSA  Web  site31),  the  ecliptic  latitude  of 
the  target  defines  when  one  can  observe  the  target  and  for  how 
long.  At  any  one  time,  Spitzer  can  observe  in  an  annulus  defined 
by  the  operational  pointing  zone,  which  can  be  conceptually 
summarized  as  “neither  too  close  nor  too  far  away  from  the 
Sun.”  It  is  about  40°  wide,  and  rotates  with  the  Sun  at  a  rate 
of  about  a  degree  a  day.  An  object  near  the  ecliptic  plane  can 
thus  only  be  observed  for  a  period  of  about  40  days  twice  a 
year;  objects  near  the  ecliptic  pole  can  be  observed  at  any  time 
during  the  year  for  as  long  as  needed.  Related  to  this,  the  IRAC 
FOV,  as  projected  onto  the  sky  for  any  given  object  on  the 
ecliptic  equator,  is  at  an  essentially  constant  angle  with  respect 
to  north  for  the  duration  of  the  ~40  day  observing  window,  and, 
~6  months  later,  is  flipped  by  180c  but  then  also  essentially 
constant  for  that  ~40  day  observing  window.  However,  the  FOV 
for  an  object  at  the  ecliptic  pole  rotates  by  about  a  degree  a  day. 
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Figure  3.  Same  as  Figure  2,  but  for  our  Orion  observations,  with  the  north-south 
boundary  extremes  separated  by  ~  1 . 6.  The  central  Chandra  (yellow)  polygon  is 
the  approximate  footprint  from  the  deep  ONC  observation  (Getman  et  al.  2005); 
the  northern  and  southern  Chandra  pointings  are  much  shallower,  and  are  taken 
from  Ramirez  et  al.  (2004b).  The  background  image  here  is  an  IRAC-2  image 
from  the  YSOVAR  campaigns  (as  opposed  to  a  cryo-era  image,  as  it  is  for  most 
of  the  other  figures  like  this). 

(A  color  version  of  this  figure  is  available  in  the  online  journal.) 

For  each  of  the  figures  showing  the  IRAC  footprints  (Fig¬ 
ures  1-13),  the  projected  outline  of  the  observation  covered  by 
the  entire  YSOVAR  data  set  is  indicated  on  top  of  a  4.5  gm 
observation  obtained  (in  most  cases)  during  the  cryogenic  era. 
The  different  regions  observed  by  the  3.6  and  4.5  gm  cameras 
are  identified;  the  nominal  target  of  the  observation  is  covered  in 
both  FOVs.  In  addition  to  the  target  of  the  observations,  there  are 
serendipitous  data  obtained  offset  from  the  target  area,  as  seen 


Figure  4.  Same  as  Figure  2,  but  for  our  Mon  R2  observations,  with  the 
north-south  boundary  extremes  separated  by  ~2 1 '.  Field  rotation,  while  present, 
is  not  as  substantial  for  this  field  as  it  is  for  others  of  our  clusters.  There  is  a 
small  amount  of  additional  X-ray  data  (not  relevant  for  this  project)  to  the  upper 
right,  and  a  portion  of  that  footprint  can  be  seen. 

(A  color  version  of  this  figure  is  available  in  the  online  journal.) 

in  the  figures.  The  footprint  of  the  serendipitously  obtained  data 
in  each  band  can  change  considerably  depending  on  the  time  of 
observation  and  the  ecliptic  latitude  of  the  target.  For  targets  at 
higher  ecliptic  latitudes,  where  field  variation  with  time  is  im¬ 
portant,  a  single  epoch  of  observation  is  indicated  in  the  figure 
in  addition  to  the  area  covered  by  the  entire  YSOVAR  data  set. 

The  approximate  ecliptic  latitude  for  each  target  is  included 
in  Table  2  as  a  rough  indication  of  how  long  the  observing 
window  was  and  how  much  field  rotation  one  can  expect  to 
see  in  these  figures.  For  Orion,  we  designed  our  observations 
to  be  less  sensitive  to  rotation;  that  plus  the  fact  that  the  Orion 
map  contains  many  FOVs  means  that  the  rotation  is  much  less 
of  a  factor  in  terms  of  how  much  of  the  sky  may  be  included 
in  both  channels  (as  opposed  to  only  one).  For  the  smaller- 
field  cluster  observations  containing  only  one  to  four  IRAC 
FOVs,  a  primary  target  field  is  monitored  in  both  IRAC  bands, 
but  substantial  serendipitous  monitoring  data  have  also  been 
obtained,  typically  in  just  one  band.  For  most  of  the  smaller- 
field  YSOVAR  cluster  targets,  such  as  NGC  1333  (Figure  2), 
the  central  region  has  complete  two-band  coverage  and  most 
of  the  objects  in  the  region  of  the  single-band  coverage  north 
and  south  of  the  target  field  have  data  over  most,  if  not  all,  of 
the  campaign.  For  L16S8  (Figure  8),  with  an  ecliptic  latitude  of 
—3°,  the  field  rotation  is  essentially  zero  during  one  ~~40  day 
window.  However,  it  was  monitored  over  more  than  one  window, 
so  the  observations  flip  by  1 80°  for  the  next  ~40  day  window. 
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Figure  5.  Same  as  Figure  2,  but  for  our  GGD  12-15  observations,  with 
the  north-south  boundary  extremes  separated  by  ~23'.  A  single  epoch  of 
observation  is  also  indicated  by  the  thinner  blue  solid  (3.6  (im)  and  red  dashed 
(4.5  /rm)  lines  to  give  an  indication  of  the  magnitude  of  the  field  rotation  for 
this  field. 

(A  color  version  of  this  figure  is  available  in  the  online  journal.) 


(The  targets  of  observation  in  Figure  8  are  the  three  regions 
with  brightest  stars  and  nebulosity  in  the  centers  of  the  three 
“stripes”  of  coverage.)  In  contrast,  for  IC  1396 A  or  Ceph  C 
(ecliptic  latitudes  of  ~60°),  the  central  5'  x  5'  region  is  covered 
throughout  the  monitoring  window,  but  the  rapid  field  rotation 
creates  a  “fan”  of  coverage  such  that  the  serendipitous  coverage 
in  the  outer  periphery  of  those  targeted  regions  have  data  in  just 
one  band,  and  only  for  a  fraction  of  the  campaign.  Some  objects 
are  monitored  in  3.6  jrm  at  the  beginning  of  the  campaign, 
and  in  4.5  jtm  at  the  end  of  the  campaign  (see  Figures  12 
and  13). 

If  one  considers  only  objects  with  light  curves  in  both 
channels  for  the  smaller-field  clusters,  even  for  the  fields  that  do 
not  rotate,  one  loses  a  significant  fraction  of  the  available  data; 
while  the  nominal  targets  of  our  observations  are  the  regions 
with  two-band  coverage,  there  are  still  good  light  curves  for 
cluster  members  in  the  regions  with  coverage  in  only  one  band. 
Generally,  for  the  smaller-field  clusters,  we  did  not  require  data 
in  both  channels,  and  instead  retained  all  light  curves  for  our 
analysis,  even  if  obtained  in  only  one  channel  (Section  3.2).  Note 
also  that  the  highest  ecliptic  latitude  among  our  targets  is  —65°; 
no  higher  ecliptic  latitude  targets  were  included,  despite  the 
possibility  of  much  longer  observing  windows,  at  least  in  part 


Figure  6.  Same  as  Figure  2,  but  for  our  NGC  2264  observations,  with  the 
north-south  boundary  extremes  separated  by  —20';  the  underlying  image  comes 
not  from  the  cryogenic  era,  but  from  the  CS1  2264  observations.  Note  that  this 
is  only  the  field  covered  as  part  of  the  original  YSOVAR  observations,  e.g., 
program  61027.  For  CSI  2264,  see  the  next  figure.  Field  rotation,  while  present, 
is  not  substantial.  The  background  image  here  is  an  IRAC-2  image  from  the 
CSI  2264  campaign  (as  opposed  to  a  cryo-era  image,  as  it  is  for  most  of  the 
other  figures  like  this). 

(A  color  version  of  this  figure  is  available  in  the  online  journal.) 


Figure  7.  Same  as  Figure  2,  but  for  our  CSI  2264  observations,  with  the 
north-south  boundary  extremes  separated  by  —0:9.  The  underlying  image 
comes  from  the  CSI  2264  observations.  The  smaller  footprints  are  the  YSOVAR- 
classic  monitored  region  from  the  previous  figure. 

(A  color  version  of  this  figure  is  available  in  the  online  journal.) 
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Figure  8.  Same  as  Figure  2,  but  for  our  L 1 688  observations,  with  the  north-south 
boundary  extremes  separated  by  ~38'.  The  field  rotation  is  essentially  non¬ 
existent  during  one  ~40  day  window,  and  then  flips  by  180  for  the  next 
~40  day  window.  The  targets  of  observation  are  the  three  regions  with  brightest 
stars  and  nebulosity  in  the  centers  of  the  three  “stripes”  of  coverage. 

(A  color  version  of  this  figure  is  available  in  the  online  journal.) 


because  we  could  not  find  any  suitable  SFRs  at  these  latitudes 
at  the  time  we  planned  our  observations.  The  field  rotation  for 
these  targets  at  ~65;  is  already  substantial. 

In  general,  because  of  the  field  rotation,  as  any  given  object 
moves  into  (or  out  from)  the  mapped  region,  the  first  (or  last) 
few  epochs  may  be  unusable  as  the  object  may  appear  in  only  a 
fraction  of  the  dithers  at  that  position  and  the  photometry  may 
thus  be  compromised  by  edge  effects.  Substantial  field  rotation 
with  time  can  have  a  significant  effect  on  faint  stars  near  bright 
stars.  Diffraction  spikes  of  bright  stars  rotate  with  time  within 
the  IRAC  FOV,  so  if  a  target  star  falls  near  a  diffraction  spike, 
this  can  introduce  a  source  of  false  longer  timescale  variability; 
see  also  Section  2.5. 

We  note  here  that  the  low  ecliptic  latitude  of  LI 688  (and 
specifically  the  180  flip  between  observing  windows)  enables 
Gunther  et  al.  (2014)  to  use  these  observations  to  characterize 
possible  artifacts  in  the  light  curve  related  to  the  position  in 
each  frame.  These  tests  confirm  that  the  statistical  criteria  for 
selecting  variable  sources  (see  Section  5  below)  is  restrictive 
enough  that  differences  between  positions  within  the  frame 
cannot  be  a  primary  contributor  for  the  sources  identified  as 
variable. 


Figure  9.  Same  as  Figure  2,  but  for  our  Serpens-Main  observations,  with 
the  north-south  boundary  extremes  separated  by  20,019,990,167  ~  23'.  Field 
rotation,  while  present,  is  not  substantial.  There  is  a  small  amount  of  additional 
X-ray  data  (not  relevant  for  this  project)  to  the  lower  right,  and  a  portion  of  that 
footprint  can  be  seen. 

(A  color  version  of  this  figure  is  available  in  the  online  journal.) 


2.4.  Warm  Spitzer  Observations:  Observing  Cadence 

More  than  95%  (524  hr)  of  our  original  YSOVAR  Spitzer 
Cycle  6  program  observing  time  was  devoted  to  a  “fast  cadence” 
mode,  designed  to  be  sensitive  to  timescales  from  ~0.15  to 
~40  days,  consistent  with  the  known  timescales  of  YSOs  due  to 
accretion-related  flickering,  rotational  modulation  of  star  spots, 
and  other  effects.  The  upper  limit  to  this  range  of  timescales  was 
set  by  the  typical  duration  of  Spitzer  visibility  windows  near  the 
ecliptic  plane,  as  discussed  above  in  Section  2.3. 

Sampling  a  light  curve  at  evenly  spaced  intervals  introduces 
period  aliasing,  and  a  bias  toward  the  detection  of  false  periods 
at  integer  fraction  multiples  of  the  sampling  interval  (see,  e.g., 
Plavchan  et  al.  2008b  or  Dawson  &  Fabrycky  2010).  This 
period  aliasing  is  common  at  integer  fraction  multiples  of  one 
day  in  ground-based  photometric  surveys  due  to  the  natural 
day-night  cycle  of  the  Earth's  rotation.  Given  typical  YSO 
rotation  periods  of  one  to  several  days,  we  therefore  chose  to 
execute  our  program  with  a  cadence  designed  to  be  compatible 
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Figure  10.  Same  as  Figure  2,  but  for  our  Serpens-South  observations,  with 
the  north-south  boundary  extremes  separated  by  ~21'.  Field  rotation,  while 
present,  is  not  substantial.  There  is  a  small  amount  of  additional  X-ray  data  (not 
relevant  for  this  project)  to  the  lower  right,  and  a  portion  of  that  footprint  can 
be  seen. 

(A  color  version  of  this  figure  is  available  in  the  online  journal.) 

with  the  Spitzer  scheduling,  and  minimize  period-aliasing  while 
retaining  sensitivity  to  variability  on  a  broad  range  of  timescales. 

The  fast  cadence  mode  for  each  region  had  a  total  time 
baseline  of  ~40  days,  varying  depending  on  the  actual  duration 
of  the  visibility  window,  with  additional  days  on  either  end 
of  the  visibility  window  for  scheduling  flexibility.  The  actual 
total  duration  of  the  fast  cadence  observations  for  each  cluster 
is  listed  in  Table  3.  In  particular,  for  AFGL  490,  Mon  R2, 
GGD  12-15,  IRAS  20050+2720,  IC  1396A,  and  Ceph  C, 
we  originally  planned  a  42  day  timespan;  for  Orion,  Serpens 
Main,  and  Serpens  South,  we  planned  a  38.5  day  timespan; 
and  for  NGC  1333,  NGC  2264  and  LI 688,  we  planned  a 
35  day  timespan.  We  divided  the  time  baseline  into  a  repeating 
set  of  3.5  day  sub-cadences  to  ease  scheduling.  Within  each 
3.5  day  period,  we  made  eight  visits,  with  time  steps  between 
visits  of  4,  6,  8,  10,  12,  14,  and  16  hr,  respectively.  By 


using  a  linearly  increasing  time  step,  we  were  able  to  evenly 
sample  in  Fourier  space  a  range  of  higher  frequency  (shorter 
timescale)  photometric  variability  than  we  would  have  been 
able  to  otherwise.  Additionally,  we  were  also  able  to  minimize 
the  total  amount  of  necessary  observing  time  to  sample  these 
high  frequencies  and  to  minimize  the  period  aliasing.  For 
comparison,  splitting  8  visits  evenly  across  a  3.5  day  cadence 
would  have  removed  our  sensitivity  to  variability  timescales 
shorter  than  ~0.5  day. 

To  further  accommodate  scheduling  flexibility,  we  also  in¬ 
cluded  a  window  of  ±2  hr  in  which  to  obtain  a  given  epoch 
of  observation  for  our  fast  cadence  observations.  This  ±2  hr 
window  effectively  resulted  in  a  randomization  of  the  cadence 
about  our  desired  times  of  observation,  which  had  the  addi¬ 
tional  benefit  of  further  reducing  period  aliasing,  in  particular 
aliases  with  periods  of  3.5  days.  We  also  allowed  for  interrup¬ 
tions  due  to  data  downlink  transfers  and  higher  priority  obser¬ 
vations  such  as  staring  mode  observations  of  transiting  exoplan¬ 
ets.  We  are  grateful  to  the  SSC  scheduling  staff  for  accommo¬ 
dating  this  fairly  complex  cadence  over  the  two  years  of  the 
survey  program. 

Overall,  we  obtained  ~40-100  epochs  per  cluster,  which  are 
visualized  in  Figure  14.  Table  2  lists  all  the  clusters,  with  in¬ 
formation  about  the  observations  (including  the  total  number  of 
epochs);  Table  3  summarizes  only  the  fast  cadence  observations, 
including  typical  values  of  the  minimum  and  maximum  size  of 
the  timestep  between  epochs,  and  the  total  time  baseline.  Our 
unevenly  spaced  epochs  enable  a  characterization  of  different 
types  of  observed  variability  on  a  variety  of  timescales;  for  fur¬ 
ther  discussion,  see  below  (Section  5.4).  While  periodic  analysis 
tools  take  advantage  of  such  unevenly  spaced  data,  Findeisen 
et  al.  (20 1 4)  point  out  that  for  auto-correlation  functions  (ACFs) 
of  non-periodic  light  curves,  timescale  sensitivity  may  be  lim¬ 
ited  by  the  largest  (not  the  smallest)  adjacent  timestep  in  the 
time  series. 

For  three  SFRs  (IC  1396A,  LI 688,  and  Ceph  C)  with  longer 
visibility  windows,  or  that  we  observed  in  multiple  visibility 
windows,  we  also  observed  in  a  “slow  cadence.”  Besides  the 
42  day  fast  cadence,  we  executed  visits  with  time  steps  between 
visits  of  1, 2,  3, 4,  5,  etc.  days  to  cover  the  remaining  duration  of 
the  visibility  window.  These  observations  constituted  ~5%  of 
our  Cycle  6  observing  program  (25  hr),  and  these  observations 
allowed  us  to  be  sensitive  to  variability  timescales  of  ~40  days 
to  ~2  yr,  e.g.,  to  objects  such  as  KH15D  (e.g.,  Winn  et  al.  2006) 
and  WL4  (Plavchan  et  al.  2008a),  as  well  as  long  term  trends  on 
timescales  of  a  year  identified  in  the  NIR  (e.g.,  Parks  et  al.  2014). 
For  LI 688  and  Ceph  C,  we  also  executed  this  slow  cadence 
observing  mode  during  the  three  other  visibility  windows 
during  our  survey.  Additionally,  Orion  had  some  slower  cadence 
follow  up,  and  NGC  2264  had  the  CSI  2264  program  at  a 
different  cadence. 

For  the  primary  statistical  analysis  discussed  in  this  paper,  we 
use  the  “standard  statistical  sample,”  e.g.,  only  the  fast  cadence 
data  (where  there  are  at  least  five  epochs  in  the  light  curve; 
see  Section  3  below).  For  about  half  of  the  clusters,  this  is  the 
entirety  of  the  data  set;  see  Tables  2  and  3,  and  the  red  points 
in  Figure  14.  A  detailed  analysis  of  the  ensemble  of  all  of  these 
observations  for  each  cluster  will  be  presented  in  papers  in 
preparation,  customized  to  each  cluster. 

For  completeness,  we  note  here  the  following  specifics  about 
those  clusters  with  data  beyond  the  YSOVAR  fast  cadence. 
LI 688,  IC  1396A,  and  Ceph  C  all  have  slow  cadence  data. 
Orion  has  slower  cadence  data  over  a  second  year.  Data  for 
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Figure  11.  Same  as  Figure  2,  but  for  our  IRAS  20050  observations,  with  the  north-south  boundary  extremes  separated  by  ~21'.  As  in  Figure  5,  a  single  epoch  of 
observation  is  also  indicated  by  the  thinner  blue  solid  (3.6  /am)  and  red  dashed  (4.5  ;rm)  lines.  Multiple  visits  of  Chandra  data  at  multiple  roll  angles  were  obtained 
here,  resulting  in  the  star-shaped  Chandra  coverage. 

(A  color  version  of  this  figure  is  available  in  the  online  journal.) 


Table  3 

Summary  of  YSOVAR  Fast  Cadence  Observations 


Cluster 

Dates  Observed 

Total  AC 
(d) 

Typical  No.  Epochs 

Typical  Min(Ar)b 
(d) 

Typical  Max(Ar)b 
(d) 

No.  Obj  with  ^5  Epochs, 
Both  Channels 

AFGL  490 

2011  Oct  19-Nov  25 

37.7 

46 

0.34 

1.95 

~600 

NGC  1333 

2011  Oct  12-Nov  14 

33.5 

72 

0.04 

1.95 

~240 

Orion"' 

2009  Oct  23-Dec  1 

39.0 

80 

0.21 

~5200 

Mon  R2 

2010  Nov  8-Dec  23 

44.4 

46 

0.26 

2.67 

-240 

GGD  12-15 

2010  Nov  16-Dec  24 

38.5 

78 

0.01 

1.82 

-340 

NGC  2264d 

2010  Nov  18-Dec  24 

35.7 

39 

0.18 

2.51 

-260 

LI  688" 

2010  Apr  12-May  17 

35.4 

80 

0.11 

0.75 

-200 

Serpens  Main 

2011  May  20-Jun  24 

34.7 

81 

0.09 

0.75 

-700 

Serpens  South 

2011  May  20-Jun  24 

34.7 

81 

0.09 

0.75 

-370 

IRAS  20050+2720 

2010  Jun  1 2-Aug  9 

57.8 

99 

0.06 

3.68 

-690 

IC  1396 A" 

2009  Sep  13-Nov  9 

56.1 

105 

0.13 

2.52 

-1480 

Ceph  Cc 

2010  Sep  18-Oct  30 

42.0 

106 

0.03 

3.68 

-400 

Notes. 

a  Total  change  in  time,  in  days,  between  first  and  last  fast  cadence  point. 

b  Typical  minimum  or  maximum  change  in  time,  in  days,  between  adjacent  points  in  the  time  series.  Also  see  Figure  26  below,  which  has  histograms 
of  the  time  step  distributions. 
c  Additional  observations  exist  in  slower  cadence. 
d  Additional  observations  are  part  of  CS1  2264. 
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Figure  12.  Same  as  Figure  2,  but  for  our  1C  1396A  observations,  with  the  north-south  boundary  extremes  separated  by  ~24'.  As  in  Figure  5,  a  single  epoch  of 
observation  is  also  indicated.  Two  pointings  of  Chandra  data  at  two  different  roll  angles  were  obtained  here,  resulting  in  the  polygon  of  Chandra  coverage. 

(A  color  version  of  this  figure  is  available  in  the  online  journal.) 


Ceph  C  over  its  fast  cadence  window  include  those  originally 
obtained  as  part  of  the  original  YSOVAR  project,  as  well  as  data 
taken  using  a  very  similar  distribution  of  time  steps  as  part  of 
a  YSOVAR-related  Chandra  program;  see  Table  2.  For  consis¬ 
tency  with  the  maximum  timescales  sampled  in  other  clusters, 
we  define  the  Ceph  C  fast  cadence  light  curve  to  be  a  window 
spanning  from  2010  September  19  through  42  days  later,  which 
is  more  representative  of  the  fast  cadence  window  in  the  other 
clusters.  We  note,  however,  that  the  sampling  within  this  window 
is  nonetheless  enhanced  relative  to  other  smaller-field  clusters, 
due  to  the  inclusion  of  extra  Spitzer  observations  obtained  to  sup¬ 
port  the  coordinated  Chandra  observations.  Two  other  clusters 
(IC  1 396A  and  IRAS  20050)  also  sample  longer  timescales  than 
the  more  typical  ~40  day  window  (see  Figure  27  below  for  an 
enlargement  of  the  fast  cadence).  We  did  not  truncate  the  last 
few  slower  cadence  points  from  the  IRAS  20050  time  series 
because  there  are  no  other  slower  cadence  observations,  and 
truncating  it  would  effectively  “orphan”  the  last  few  points.  The 
IC  1396A  cadence  does  not  provide  a  clean  breakpoint,  so  it 
also  has  points  covering  a  slightly  longer  window  than  “typi¬ 
cal”  included  in  its  fast  cadence  (see  Table  3  for  the  total  time 
in  the  fast  cadence  for  all  clusters).  IC  1396A  was  one  of  the 
first  SFRs  to  be  monitored  intensively  with  Spitzer  (Morales- 
Calderon  et  al.  2009),  so  it  has  additional  cryo-era  monitoring 


not  shown  in  Figure  14.  Serpens  Main  also  has  some  cryo-era 
monitoring.  GGD  12-15  and  L1688  also  have  staring  data,  not 
included  here. 

2.5.  Warm  Spitzer  Observations:  Data  Reduction 

We  started  with  the  IDL  package  Cluster  Grinder  (G09), 
which  has  been  used  by  several  other  projects  (e.g.,  Megeath 
et  al.  2012),  and  then  made  custom  modifications  to  make  it 
suitable  for  this  time  series  data  set.  We  started  with  the  basic 
calibrated  data  (BCD)  images  released  by  the  SSC.  The  BCD 
images  used  correspond  largely  to  software  version  S 1 8. 1 8  (with 
about  20%  S 1 9.0  and  S 1 9. 1 );  earlier  reductions  such  as  those  in 
MCI  1  used  a  mix  of  SI  8. 12,  18.14,  and  18.18.  Each  BCD  frame 
was  processed  for  standard  bright  source  artifacts  and  combined 
into  mosaics  (at  CK'86  per  pixel  grid  resolution)  by  exposure  time 
(long  or  short),  epoch,  and  bandpass.  Cosmic  ray  hits  and  other 
transient  artifacts  are  flagged  during  mosaic  construction  using 
redundant  data  at  each  pixel  position  in  the  final  mosaic  grid  as  a 
reference  for  the  nominal  value  and  allowable  internal  variation. 

Since  the  HDR  mode  obtains  a  long  and  a  short  frame  at 
each  pointing,  we  need  to  combine  data  from  the  two  expo¬ 
sures.  HDR  mosaics  are  merged  together  on  a  pixel-by-pixel 
basis,  with  appropriately  scaled  short-frame  HDR  values  above 
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Figure  13.  Same  as  Figure  2,  but  for  our  Ceph  C  observations,  with  the  north-south  boundary  extremes  separated  by  ~20'.  As  in  Figure  5,  a  single  epoch  of  observation 
is  also  indicated.  The  field  rotation  here  is  quite  significant. 

(A  color  version  of  this  figure  is  available  in  the  online  journal.) 


2008  2010  2012  2010.0  2010.5  2011.0  2011.5  2012.0 


Figure  14.  Visual  overview  of  the  observing  seasons  and  the  cadence  for  all  clusters  in  the  YSOVAR  project.  On  the  left  side,  all  of  the  YSOVAR  data  are  depicted 
(where  Spitzer  s  cryogen  ran  out  2009  May  15,  MJD  54966.925).  On  the  right  side  is  an  expanded  view  of  the  data  highlighted  here  and  used  for  the  “standard 
set  for  statistics”  (see  Section  3.2).  The  designated  “fast  cadence”  observations  are  shown  in  red,  with  other  YSOVAR-classic  observations  shown  in  black.  Other 
YSOVAR-affiliated  programs  are  shown  in  blue;  the  large  post-cryo  monitoring  of  NGC  2264  is  CSI  2264,  and  the  cryogenic  observations  of  IC  1396A  and  Serpens 
are  also  shown.  Each  individual  observation  is  marked  by  a  “|”  symbol  in  which  the  length  of  the  line  is  proportional  to  the  number  of  observations  in  a  window  12  hr 
before  and  after  that  particular  observation.  When  the  observations  are  so  dense  that  the  symbols  overlap,  the  frequency  of  observations  can  thus  be  judged  from  the 
thickness  of  the  line.  An  enlargement  of  the  red,  fast-cadence  sequences  can  be  found  in  Figure  27  below. 

(A  color  version  of  this  figure  is  available  in  the  online  journal.) 
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a  set  threshold  (individual  pixel  values  corresponding  to  sources 
MOth  mag,  dependent  on  background  level)  replacing  signifi¬ 
cantly  compromised  pixels  in  the  long-frame  HDR  mosaic. 

Point  source  detection  and  aperture  photometry  are  then 
performed  for  each  channel  and  each  AOR/epoch.  The  aperture 
radius  was  2"4,  with  a  sky  annulus  having  inner  and  outer  radii  of 
2'.'4  and  7"2,  respectively.  We  adopted  Vega  standard  magnitude 
zero  points  of  19.30  and  18.64  for  1  DN  s^1  total  flux  at  3.6  and 
4.5  /ira,  respectively.  These  values  include  standard  corrections 
for  our  chosen  aperture  and  sky  annulus  sizes  (Reach  et  al.  2005; 
Carey  et  al.  2010). 

Extensive  tests  were  run  on  aperture  photometry  obtained 
in  three  different  ways.  Our  initial  approach  was  to  obtain 
photometry  from  the  mosaic  created  for  each  AOR,  with  the 
assigned  time  for  that  point  being  the  average  time  for  the  AOR. 
We  also  obtained  photometry  on  each  BCD  frame  individually, 
with  the  individual,  specific  time  of  that  observation  being 
assigned  to  that  point.  Finally,  we  took  the  set  of  all  individual 
BCD  measurements  for  a  given  source  within  a  given  AOR, 
and  took  the  mean  brightness  and  the  mean  time  for  that  set 
of  observations.  This  latter  approach  proved  to  have  the  best 
quality  photometry.  It  facilitated  the  treatment  of  two  well- 
known  systematic  effects  (residual  gain  and  pixel  phase  effects) 
while  keeping  our  signal-to-noise  ratios  (S/Ns)  close  to  ideal.  It 
provided  the  highest  S /N  light  curves,  even  if  the  time  resolution 
is  slightly  lower  than  what  might  be  theoretically  possible. 
Effectively,  this  choice  means  that  a  S/N~5  is  required  for 
a  viable  source  in  a  given  BCD,  such  that  the  net  S/N  of 
the  source  over  the  entire  AOR  is  at  least  ~10.  If  a  source 
fell  on  the  edge  of  the  BCD  such  that  the  two  native  pixel 
radius  was  not  entirely  within  a  given  BCD,  no  measurement 
is  obtained.  Measurements  were  also  excluded  as  non-viable  if 
the  measurement  (in  the  aperture  or  annulus)  included  a  masked 
IRAC  pixel  such  that  obtaining  a  finite  measurement  was  not 
possible,  or  if  the  source  was  faint  enough  (or  the  sky  variation 
high  enough)  that  the  net  measured  flux  was  <0.  If  there  were 
three  or  more  valid  detections  of  a  given  source  in  a  given  AOR, 
we  could  perform  outlier  rejection,  and  outliers  were  flagged; 
the  rest  of  the  measurements  were  combined  by  uniform  weight 
arithmetic  mean.  The  uncertainties  were  added  in  quadrature 
and  divided  by  the  included  measurement  count. 

Source  matching  between  epochs  and  wavelengths  was  per¬ 
formed  by  position,  taking  the  nearest  source  within  1".  Final 
positions  for  each  source  were  the  mean  of  the  individual  mea¬ 
surements.  Our  astrometric  residuals  as  compared  to  2MASS 
suggest  an  uncertainty  of  <200  mas  root-mean-square  (rms), 
after  a  single  iteration  refinement  of  the  astrometry  using  a 
first  set  of  mosaics  constructed  blindly  with  the  BCD-delivered 
world  coordinate  system  (WCS).  Objects  with  fewer  than  five 
viable  detections  over  separate  epochs  in  a  single  band  were  not 
retained  as  valid  light  curves. 

We  note  here  that  there  are  some  residual  instrumental 
column  pulldown  effects  (see  the  IRAC  Instrument  Handbook 
for  more  information)  in  some  areas  near  bright  sources.  This  is 
a  more  significant  problem  for  clusters  with  a  large  number 
of  very  bright  sources  (e.g.,  Mon  R2;  see  brightnesses  in 
Figures  1-13).  The  locations  of  these  artifacts  can  change  with 
time,  particularly  in  fields  at  high  ecliptic  latitude  and  thus 
with  significant  field  rotation  during  a  given  visibility  window 
(see  Section  2.3),  and  can  have  particularly  significant  effects 
on  sources  fainter  than  —  1 2th  mag.  Objects  falling  in  such 
regions  that  are  faint  and  that  had  non-repeating  structure  in  their 
light  curves  should  be  particularly  scrutinized  (including  visual 


inspection  of  the  input  frames)  for  the  validity  of  their  light 
curves.  Points  obviously  compromised  by  these  effects  should 
additionally  be  identified  as  non-viable  points  and  rejected. 
Discussions  of  this  process  for  individual  objects  will  appear 
in  the  clusters  papers  to  follow. 

Even  after  all  of  the  processing  to  this  point,  some  occasional 
outliers  in  the  light  curves  remain.  To  identify  outliers  within  the 
fast  cadence  data,  we  initially  identified  points  several  a  away 
from  the  mean  for  each  light  curve.  We  used  this  approach  in 
MCI  1,  with  a  different  (much  earlier)  processing  of  the  original 
Spitzer  data.  However,  our  more  recent,  improved  processing 
described  here  reduced  the  number  of  outliers.  We  experimented 
with  algorithms  for  rejecting  outliers  from  each  light  curve, 
but  were  unable  to  identify  a  technique  that  rejected  visually 
spurious  points  while  preserving  apparent  bona  fide  variability. 
An  approach  as  we  used  in  MCI  1  would,  for  example,  exclude 
points  at  the  beginning  and  end  of  a  light  curve  with  a  long 
trend.  Our  estimates  suggest  that  implementing  an  outlier 
rejection  strategy  affects  <10%  of  the  light  curves  tagged  as 
intrinsically  variable  (as  per  methods  described  further  below), 
so  we  retained  all  points  in  the  light  curves.  Outliers  may  still  be 
identified  in  a  few  specific  individual  cases;  those  outliers  are 
then  effectively  excluded  from  the  light  curve,  and  are  discussed 
when  relevant. 

2.6.  Warm  Spitzer  Observations:  Noise  Floor 

Many  tests  of  variability  (e.g.,  a  x2  test;  see  Section  5.2) 
depend  on  accurate  error  estimates  for  each  data  point  in  the 
time  series.  Initial  estimates  of  uncertainties  were  derived  by 
combining  three  terms  in  quadrature  (see  also  the  discussion  in 
Megeath  et  al.  2004);  shot  noise  in  the  aperture,  shot  noise  in 
the  mean  background  flux  per  pixel  integrated  over  the  aperture, 
and  the  standard  deviation  of  the  sky  annulus  pixels  to  account 
for  the  influence  of  non-uniform  nebulous  background.  Such 
error  values  proved  to  be  a  slight  underestimate  of  the  true 
uncertainty,  and  do  not  take  into  account  a  “floor”  in  the  errors 
from  calibration  errors,  primarily  the  intrapixel  gain  variation. 

We  wanted  to  find  a  value  for  this  floor  that  worked  for 
all  clusters  (since  all  of  the  mapping  observations  discussed 
here  were  obtained  in  the  same  way).  We  also  wanted  to  be 
conservative  in  that  we  wanted  the  set  of  objects  we  selected 
as  variable  at  the  end  of  this  process  to  be  reliably  variable 
(with  very  few,  if  any,  non-variables  selected  as  variable),  even 
if  it  meant  that  our  sample  would  not  be  complete  in  that 
some  legitimately  (lower  level)  variable  objects  would  not  be 
selected.  In  other  words,  we  wished  to  err  slightly  on  the  side 
of  overestimating  errors.  We  investigated  the  rms  error  for  each 
light  curve  for  each  object  versus  magnitude  and  versus  the 
median  photometric  uncertainty  for  all  sources.  In  MCI  1,  we 
estimated  the  true  error  as  a  function  of  magnitude  to  be  single¬ 
valued  as  a  function  of  magnitude  bin.  In  general,  this  means  that 
the  errors  will  be  correct  on  average,  but  there  will  be  specific 
objects  (e.g.,  those  in  high  background  regions)  for  which  they 
are  too  low. 

Here,  we  aimed  to  improve  on  this  estimate  by  finding  a 
value  for  the  noise  floor  that  can  be  added  in  quadrature  to  the 
error  estimates  for  each  point,  preserving  the  errors  obtained 
for  each  point,  localized  in  time  and  space.  For  each  channel, 
for  every  object  in  each  cluster  having  at  least  20  epochs  in  its 
light  curve  so  as  to  ensure  well-determined  values,  we  plotted 
the  measured  rms  (cr)  of  the  light  curve  against  its  median 
photometric  uncertainty;  see  Figure  15.  Orion  alone  represents 
about  one-quarter  of  the  light  curves  shown  in  Figure  15,  so  it 
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Figure  15.  Plots  of  measured  rms  (cr)  vs.  median  photometric  uncertainty  for  1 1  (left)  and  12  (right).  The  top  row  is  for  all  clusters  except  Orion,  the  middle  row  is 
Ceph  C  only,  and  the  bottom  row  is  Orion  only.  The  straight  red  line  is  the  one-to-one  relationship  between  these  parameters,  e.g.,  the  expected  relationship  if  there 
was  no  “noise  floor";  the  curved  red  line  is  the  empirically  derived  curve  which  we  used  to  determine  a  noise  “floor"  of  0.01  mag  for  II  and  0.007  mag  for  12,  which 
we  then  added  in  quadrature  to  the  individual  errors  obtained  for  each  point.  To  appear  in  this  plot,  objects  must  have  more  than  20  epochs.  Ceph  C  appears  separately 
to  give  an  indication  of  what  these  plots  look  like  for  individual  clusters.  We  combined  all  the  clusters  together  to  better  determine  the  empirical  floor,  even  though 
there  are  some  variations  among  the  clusters;  see  the  text. 

(A  color  version  of  this  figure  is  available  in  the  online  journal.) 


is  useful  in  some  contexts  to  separate  Orion  from  the  analysis 
of  the  rest  of  the  ensemble.  Data  for  an  individual  smaller-field 
cluster  are  thus  shown  in  the  figure  for  reference. 

The  largest-amplitude  sources  in  Figure  15  are  true  variables. 
To  determine  the  noise  floor,  we  are  interested  in  how  the 
observed  rms  for  non-variables  over  many  epochs  compares 
to  the  formal  errors.  The  bulk  of  the  distribution  of  points 
in  Figure  15,  e.g.,  the  non-variable  and  less  noisy  sources  in 
the  ensemble  of  points,  fall  in  the  same  region  in  each  plot; 
the  distributions  run  smoothly  down  to  an  asymptotic  value. 
However,  there  are  three  components  contributing  to  the  overall 
scatter  of  the  regions  with  more  objects  seen  in  Figure  15, 
where  they  are  not  necessarily  seen  as  distinct  features.  As  we 
move  to  brighter  objects  (generally  those  with  smaller  median 
photometric  uncertainty),  we  have  a  statistically  larger  chance 
that  the  objects  will  be  cluster  members,  e.g.,  legitimately  young 
and  therefore  expected  to  have  a  greater  likelihood  of  being 
highly  variable.  There  is  a  transition  between  the  short  and  long 


HDR  frames  (at  ~10th  mag)  which  affects  the  measured  error 
for  sources  that  are  faint  in  the  short  frames  but  not  in  the  long 
frames.  There  are  variations  in  depth  (~  1 6th— 1 7th  mag)  reached 
across  clusters  due  to  variable  and  nebulous  backgrounds  (the 
variations  in  depth  can  also  be  seen  in  the  faint  limits  of 
Figures  19  and  20,  discussed  below)  such  that  a  given  median 
photometric  uncertainty  need  not  necessarily  refer  to  stars  of 
roughly  the  same  brightness  among  the  clusters.  Figure  15 
includes  an  individual  cluster  (Ceph  C)  to  give  an  indication  of 
this  cluster-dependent  variation  in  the  floor;  it  can  be  seen  that 
fewer  points  fall  below  the  red  line  on  the  left  side  of  each  panel 
in  Ceph  C  than  they  do  in  the  Orion  plots  on  the  bottom  row. 

Close  inspection  of  Figure  1 5  reveals  that  the  asymptotic  noise 
floor  value  is  on  average  slightly  higher  for  the  3.6  /im  channel 
and  slightly  lower  for  the  4.5  jrm  channel.  This  is  consistent  with 
expectations  that  the  dominant  source  of  error  in  the  asymptotic 
noise  floor  is  from  residual  intrapixel  gain  variations,  since  the 
intrapixel  gain  effect  is  lower  in  12  than  in  II.  There  are  also 
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Figure  16.  Plots  of  measured  rms  (cr)  vs.  mean  magnitude  for  11  (left)  and  12  (right).  The  top  row  is  for  all  clusters  except  Orion,  the  middle  row  is  L1688  only,  and 
the  bottom  row  is  Orion  only.  There  is  more,  and  more  significant,  variation  among  the  brighter  sources.  (Some  very  large  rms  values  fall  above  the  range  shown  in 
this  plot;  the  range  is  limited  here  for  clarity.) 


more  points  with  significantly  larger  a  in  the  4.5  yum  channel. 
The  4.5  /rm  channel  has  more  contributed  noise  in  any  given 
light  curve — shocks  and  scattered  light,  which  are  very  common 
in  our  target  regions,  contribute  to  the  overall  scatter  in  a 
at  4.5  /iin.  There  is  also  some  lower-level  contribution  from 
polycyclic  aromatic  hydrocarbon  (PAH)  features  (also  common 
in  these  regions)  in  3.6  jtm  (see,  e.g.,  Flagey  et  al.  2006).  These 
contributions  to  both  channels  mean  that  there  is  some  cluster-to- 
cluster  variation  in  the  uncertainties  associated  with  individual 
sources  at  a  given  magnitude. 

We  determined  that  a  systematic  error  of  0.01  mag  in  II 
and  0.007  mag  in  12  (both  indicated  in  Figure  15)  provides 
the  best  fit  to  the  photometric  rms  seen  in  each  channel  in 
that  it  is  a  conservative  representation  of  the  distribution  as 
the  distributions  approach  the  asymptote.  These  values  are  also 
quite  comparable  to  the  magnitude  of  the  intrapixel  gain  effect 
in  each  of  these  channels,  and  it  is  not  expected  that  the  errors  in 
IRAC  photometry  from  mapping  observations  would  be  less 
than  0.1%.  Very  similar  values  are  obtained  via  a  slightly 
different  approach  for  the  CSI  2264  mapping  in  Cody  et  al. 
(2014),  where  the  populations  of  members  and  non-members  are 


better  understood.  We  added  our  empirically  derived  systematic 
values  in  quadrature  to  each  of  the  individual  errors  obtained  by 
our  pipeline. 

Figure  16  plots  the  light  curve  rms  (ct)  derived  from  the  cor¬ 
rected  light  curves  against  the  mean  magnitude  for  that  light 
curve.  Here,  again,  one  can  see  more,  and  more  significant, 
variation  among  the  brighter  sources.  Relatively  few  sources 
have  much  larger  uncertainties  because  they  are  close  to  bright 
neighbors,  chip  edges,  or  affected  by  other  instrumental  ar¬ 
tifacts;  many  of  the  sources  with  large  rms  are  legitimately 
variable.  For  sources  brighter  than  13th  mag,  the  total  uncer¬ 
tainty  is  dominated  by  the  error  floor.  The  overall  noise  in¬ 
creases  significantly  fainter  than  about  14th  magnitude,  and 
the  rapid  falloff  of  the  number  of  objects  detected  in  the  sur¬ 
vey  beyond  16th  magnitude  is  immediately  apparent  (see  also 
Figures  19  and  20  below).  There  are  small  variations  between 
clusters,  which  are  caused  by  different  levels  of  diffuse  emis¬ 
sion  and  a  different  degree  of  instrumental  artifacts  from,  e.g., 
residual  pull-down.  A  sample  individual  cluster  is  included  in 
Figure  16  as  an  indication  of  what  is  seen  in  the  individual 
smaller-field  clusters. 
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Table  4 

Summary  of  Cryo-epoch  Spitzer  Observations 


Cluster 

AORKEYa 

Obs.  Date 

Adop.  Date  b 

A(f)c 

Notes 

AFGL  490 

3654656,3663104 

2004  Oct  8,  2004  Sep  24 

2004.9 

6.9 

Reduction  from  G09 

NGC  1333 

3652864,5793280,  4316672 

2004  Feb  10,  2004  Sep  8,  2004  Feb  3 

2004.5 

7.3 

Reduction  from  G09,  Gutermuth 
et  al.  (2008b);  some  data  from 

Orion 

(Many) 

2004  Mar-2004  Oct 

2004.5 

5.9 

c2d  included  where  necessary 
Reduction  from  Megeath  et  al. 

Mon  R2 

3659008,  3668224 

2004  Mar  13,  2005  Apr  7 

2004.3 

6.6 

(2012) 

Reduction  from  GOO 

GGD  12-15 

3658752,  3667968 

2004  Mar  31,  2004  Mar  15 

2004.3 

6.6 

Reduction  from  GOO 

NGC  2264 

3956480,  3956736,  3956992, 

2004  Mar  6,  2004  Oct  8,  2004  Mar  6, 

2004.5 

6.4 

Re-reduced  in  G09  style 

3957248 

2004  Oct  8 

L1688 

3652096,  4321280 

2004  Mar  7,  2004  Aug  26 

2004.3 

6.7 

Reduction  from  G09;  some  data 
from  c2d  included  where 

Serpens  Main 

3652352,3661568,  12649216 

2004  Apr  1 , 2004  Apr  6,  2005  Apr  1 2 

2004.3 

7.2 

necessary 

Reduction  from  G09;  some  data 
from  c2d  included  where 

Serpens  South 

19957248,  19958016,  19998464, 

2006  Oct  27,  2006  Oct  27,  2006  Oct 

2006.9 

4.6 

necessary 

Reduction  from  Gutermuth  et  al. 

19998720,  20002304,  20002560, 

28,  2006  Oct  28,  2006  Oct  27,  2006 

(2008b) 

20018688,20018944 

Oct  27,  2006  Oct  28,  2006  Oct  28 

IRAS 

20050+2720 

3656448,  3665152 

2004  Jun  9,  2004  Sep  25 

2004.5 

6.1 

Reduction  from  G09 

IC  1396A 

3959040,  4316416 

2003  Dec  20,  2004  Jun  23 

2004.0 

6.4 

Re-reduced  in  G09  style 

Ceph  C 

3655936,  3664384 

2003  Dec  19, 2003  Dec  19 

2004.0 

6.4 

Reduction  from  G09 

Notes. 

a  An  AOR  is  an  Astronomical  Observation  Request,  the  fundamental  unit  of  Spitzer  observing.  An  AORKEY  is  the  unique  eight-digit  integer  identifier  for  the  AOR, 
which  can  be  used  to  retrieve  these  specific  data  from  the  Spitzer  Heritage  Archive. 

b  This  is  the  adopted  date  of  the  IRAC  observations  used  in  the  present  paper  for  the  cryo  epoch  (e.g.,  the  time  that  was  used  to  represent  the  date  of  those  observations 
taken  over  more  than  one  epoch).  Note  that  two  clusters  have  monitoring  observations  taken  during  the  cryo  era,  but  this  —single-epoch  earliest  cryo  observation  is 
what  is  used  in  the  present  paper  as  the  “cryo  epoch”  for  comparison  to  the  “YSOVAR  epoch." 
c  A((),  in  years,  between  adopted  cryo  epoch  and  mean  YSOVAR  fast-cadence  time. 


2. 7.  Cryogenic-era  Spitzer  Data 

Early  in  the  Spitzer  mission,  each  of  our  12  regions  was 
observed,  often  as  part  of  the  guaranteed  time  observations  or 
the  original  Cores-to-Disks  (c2d)  Legacy  program  (Evans  et  al. 
2003,  2009b);  these  observations  are  summarized  in  Table  4. 
In  most  cases,  the  data  were  obtained  at  two  epochs.  For 
clusters  near  the  ecliptic  plane,  asteroids  that  happen  to  be 
passing  through  the  region  at  the  time  of  observation  appear 
as  long  wavelength  sources  that  can  resemble  embedded  low- 
mass  YSOs.  Thus,  the  two  epochs  were  often  planned  to  be 
separated  by  a  time  of  the  order  of  hours  (e.g.,  the  earliest  epochs 
of  Serpens,  much  of  LI 688)  to  allow  the  moving  solar  system 
objects  enough  time  to  move,  and  the  corresponding  pixels  to  be 
rejected  as  outliers  when  the  individual  frames  were  combined. 
For  some  observations  (e.g.,  Orion),  to  facilitate  bright  source 
artifact  mitigation,  the  epochs  were  separated  by  ~6  months 
such  that  the  orientation  of  the  arrays  in  the  later  observation 
is  180°  from  that  in  the  early  observation  (see  Section  2.3 
above).  For  the  most  robust  detections,  the  cryogenic  data  can  be 
combined  into  a  single  early  epoch,  or  one  can  compare  the  two 
epochs  to  constrain  variability  on  timescales  of  months  (e.g., 
Megeath  et  al.  2012  for  Orion  at  IRAC  bands)  or  hours  (e.g.. 
Rebull  et  al.  2007  for  Perseus  at  24  pm). 

To  retain  information  from  the  individual  epochs,  the  cryo¬ 
genic  data  for  all  of  our  clusters  were  reduced  identically  to  the 
YSOVAR  monitoring  data,  except  using  cryogenic  calibrations 
(e.g.,  the  zero  point  and  threshold  between  the  short  and  long 
HDR  frames).  Source  detection  was  performed  independently 
in  each  band  for  each  epoch.  However,  for  reasons  we  now  de¬ 


scribe,  measurements  derived  from  individual  epochs  from  the 
cryogenic  era  cannot  be  used  generally  over  the  entire  region, 
and  the  identical  approach  adopted  for  the  rest  of  the  YSOVAR 
data  cannot  be  used  generally  for  the  cyrogenic  data. 

For  the  YSOVAR  data,  each  epoch  covered  the  region  of 
interest  with  sufficient  redundancy,  so  that  data  reduction  can 
be  performed  on  a  per-epoch  basis  as  described  in  Section  2.5 
above.  For  some  of  the  cryogenic  era  data,  the  observations 
were  designed  to  have  sufficient  redundancy  in  the  region  of 
interest  only  once  all  of  the  cryogenic  era  observations  were 
combined.  Specifically,  for  some  of  those  clusters  that  were 
observed  in  two  epochs,  for  moving  object  identification,  the 
observations  from  a  single  epoch  were  not  designed  to  robustly 
measure  objects  at  that  epoch;  it  was  envisioned  that  the  two 
epochs  would  be  combined  to  create  a  net  mosaic  of  the  inertial 
targets.  For  example,  one  epoch  might  have  two  frames  with 
another  two  frames  at  a  later  epoch.  Combining  those  after  the 
fact  provides  four  frames  of  dithered  observations  to  remove 
artifacts  and  moving  objects.  However,  a  single  epoch,  having 
only  two  frames,  does  not  necessarily  have  enough  data  for 
robust  photometry.  A  complication  comes  from  the  combination 
of  individual  frames  to  compose  a  map  at  a  given  epoch.  The 
maps  at  any  one  epoch  are  constructed  from  multiple  pointings, 
with  overlap  between  each  pointing.  Therefore,  in  those  thin 
regions  between  adjacent  pointings,  there  are  more  than  two 
frames  at  a  given  epoch.  Because  our  pipeline  requires  at  least 
three  frames  per  position,  it  only  derives  (and  archives)  high- 
quality  photometry  at  either  of  those  two  epochs  for  that  thin 
region  between  adjacent  pointings  in  the  map,  not  the  entire  map 
at  that  epoch.  Therefore,  if  one  downloads  the  earliest  epoch  data 
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from  our  archive  for  those  clusters,  one  obtains  a  “grid-shaped 
web”  of  viable  data,  not  complete  coverage  of  the  region.  We 
have  retained  all  of  these  individual  epoch  data  in  the  data  we 
deliver  to  IRSA  because  it  may  be  useful  to  future  researchers  to 
have  individual  measurements  of  sources  serendipitously  falling 
in  those  regions.  However,  for  users  not  aware  of  this  “feature,” 
it  will  seem  strange  to  only  have  cryogenic  data  over  a  web  of 
coverage  in  those  early  epochs. 

In  the  context  of  the  present  paper,  therefore,  all  of  the  earliest 
cryogenic-era  Spitzer  observations  were  combined  into  one 
mosaic  per  epoch  per  channel,  which  was  then  run  through 
Cluster  Grinder.  The  time  assigned  to  this  early  epoch  of 
observation  is  the  average  time  of  all  of  the  observations  that 
went  into  the  mosaic,  summarized  in  Table  4. 

As  previously  stated,  IC  1396 A  and  Serpens  Main  were 
explicitly  monitored  in  the  cryogenic  era  with  Spitzer,  and 
these  data  were  also  reprocessed,  retaining  individual  epochs 
of  observation. 

The  data  were  bandmerged  across  Spitzer  bands  by  position, 
within  a  search  radius  of  1".  (Recall  that  the  IRAC  photometry 
apertures  we  used  were  2'.'4.) 

Gutermuth  et  al.  (2008b,  2009,  2010)  present  methodology 
for  identifying  YSOs  from  the  cryogenic  catalog.  The  details  of 
the  selection  process  appear  in  those  papers,  but  in  summary, 
multiple  cuts  in  multiple  color-color  and  color-magnitude 
diagrams  are  used  to  identify  YSO  candidates,  as  distinct 
from,  e.g.,  extragalactic  and  nebular  contamination.  This  color 
selection  is  part  of  Cluster  Grinder,  and  thus  we  have  this 
classification,  based  on  the  cryogenic  data,  for  all  of  our  clusters; 
we  used  it  in  Section  2.1.2  and  in  part  of  Table  1.  We  have 
adopted  this  YSO  selection  mechanism  as  part  of  one  of  the 
primary  YSOVAR  sample  definitions;  see  Section  3. 

2.8.  New  Data  at  Optical  and  NIR  Wavelengths 

For  completeness,  we  note  here  that,  in  addition  to  the  new 
IRAC  data,  we  often  also  obtained  contemporaneous  optical 
and  NIR  observations  from  the  ground  using  a  variety  of 
telescopes.  The  details  of  each  telescope/camera  and  set  of 
accompanying  observations  is  beyond  the  scope  of  this  paper, 
and  is  (or  will  be)  provided  in  each  paper  discussing  the  results 
from  each  cluster.  For  example,  MCI  1  discussed  observations 
obtained  from  four  other  ground-based  telescopes  obtained  in 
2009-2010,  contemporaneous  with  the  first  epoch  of  Orion 
observations;  Cody  et  al.  (2014)  discussed  observations  obtained 
contemporaneously  with  the  CSI  2264  campaign.  Because  these 
ancillary  monitoring  data  are  different  for  each  cluster,  we  omit 
them  here  for  clarity. 

2.9.  Other  IR  Archival  Data 

We  also  included  2MASS  JHKS  data  for  all  of  our  clusters.  For 
two  clusters  (NGC  1333,  LI 688),  deeper  than  survey  2MASS 
data  (the  6x  data)  are  available  from  the  2MASS  archive, 
and  are  included  where  available.  Such  detections  were  also 
bandmerged  to  the  rest  of  the  catalog  by  position,  within  a 
search  radius  of  1". 

For  three  of  our  clusters  (LI 688,  Mon  R2,  Serpens  Main), 
deeper  NIR  data  from  the  UKIRT  Infrared  Deep  Sky  Survey 
(UKIDSS;  Lawrence  et  al.  2007)  broadband  data  were  publicly 
available,  but  not  necessarily  at  all  or  even  most  of  the  UKIDSS 
bands  (Z,  Y,  J,H,K)\  in  all  three  cases,  there  were  at  least  X-band 
data,  which  is  important  for  our  calculation  of  SED  class  (see 
Appendix  B).  Where  available,  these  data  are  also  bandmerged 


to  the  rest  of  the  catalog  by  position,  within  a  search  radius 
of  1".  Details  of  which  bands  are  available  and  to  what  depth 
will  be  included  in  the  individual  papers  associated  with  the 
relevant  clusters. 

For  two  of  our  clusters,  AFGL  490  and  Ceph  C,  data  in  two 
Spitzer  bands  are  available  from  one  of  the  warm  portions  of  the 
Spitzer  program  called  the  Galactic  Legacy  Infrared  Mid-Plane 
Survey  Extraordinaire  (Benjamin  et  al.  2003).  Where  available, 
those  single-epoch,  shallow  data  were  also  bandmerged  to  the 
rest  of  the  catalog  by  position,  within  a  search  radius  of  1".  Since 
those  measurements  are  independent  measures  of  these  objects 
at  [3.6]  and  [4.5],  these  measurements  were  retained  as  distinct 
points  from  the  cryo  [3.6]  and  [4.5],  or  the  means  of  our  light 
curves. 

For  three  other  clusters,  NGC  1333,  LI 688,  and  Serpens 
Main,  data  are  available  from  the  c2d  (Evans  et  al.  2003, 
2009b)  program  data  deliveries.  The  data  used  for  these  de¬ 
liveries  are  typically  the  same  BCDs  as  were  used  in  the  cryo¬ 
genic  data  (Section  2.7).  As  such,  then,  they  are  not  indepen¬ 
dent  measurements,  and  these  data  were  only  used  to  supple¬ 
ment  our  cryogenic-era  catalog  if  a  band  was  missing,  e.g., 
because  G09  had  identified  it  as  having  insufficient  S/N  in  that 
reduction. 

Naturally,  each  cluster  has  different  amounts  of  additional 
data  in  the  literature,  and  the  details  of  exactly  which  data  are 
included  (and  how  it  was  merged  to  the  rest  of  the  catalog)  will 
appear  in  the  corresponding  YSOVAR  cluster  paper. 

2.10.  Chandra  Data 

X-ray  observations  are  an  excellent  complement  to  mid-  and 
near-IR  observations  in  regions  of  star  formation.  Due  to  the 
strong  correlation  between  X-ray  luminosity  and  age.  X-ray 
emission  is  particularly  effective  in  identifying  YSOs  that  are 
free  of  IR  excess  (Class  Ills),  such  that  X-ray  surveys  provide 
samples  of  young  stars  unbiased  by  their  IR  characteristics. 

X-ray  observations  can  penetrate  up  to  Ay  ~  500  mag  into 
a  star  forming  cloud  with  very  deep  integrations  (Grosso  et  al. 
2005).  However,  even  with  shallow  integrations,  one  can  reach 
as  deep  or  deeper  in  the  X-rays  than  many  NIR  surveys  in 
the  JHK  bands.  The  Chandra  X-ray  Observatory,  with  its 
high  angular  resolution  mirrors  and  low-noise  detectors,  is 
particularly  effective  in  resolving  crowded  fields  down  to  0!5 
scales.  Furthermore,  in  X-rays,  OB  stars  are  often  not  much 
brighter  than  pre-main-sequence  stars,  so  close  companions, 
even  when  associated  with  OB  stars,  can  be  identified  (Stelzer 
et  al.  2005).  Chandra  data  (or  any  other  X-ray  data,  e.g., 
from  the  European  Space  Agency’s  XMM-Newton  observatory) 
and  Spitzer  data  combined  have  been  shown  to  yield  a  more 
complete  survey  of  members;  there  are  many  examples  of  this 
in  the  literature.  For  some  of  the  clusters  in  our  target  list,  an 
incomplete  listing  of  such  work  could  include  Winston  et  al. 
(2007,  2010)  for  NGC  1333  and  Serpens  Main,  respectively, 
Gunther  et  al.  (2012)  for  IRAS  20050,  Getman  et  al.  (2006)  for 
IC  1396A,  or  Imanishi  et  al.  (2001)  for  L1688. 

Identification  of  additional  cluster  members  is  the  primary 
purpose  for  which  we  include  X-ray  data,  where  available,  for 
the  YSOVAR  clusters.  We  can  then  compare  variability  char¬ 
acteristics  for  the  X-ray-detected  sample  with  that  from,  e.g., 
the  IR-selected  sample.  With  the  exception  of  AFGL  490,32  all 
of  the  YSOVAR  clusters  have  been  observed  by  Chandra  us¬ 
ing  the  Advanced  CCD  Imaging  Spectrometer  for  wide-field 


32  As  of  early  2014,  AFGL  490  has  not  been  observed  by  XMM-Newton  either. 
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imaging  (ACIS-I),  and  nearly  all  of  them  have  already  been 
published.  The  date,  duration,  and  aimpoints  of  these  observa¬ 
tions  are  listed  in  Table  5,  along  with  citations  to  the  literature 
where  possible;  footprints  of  the  observations  are  included  in 
Figures  2-13. 

In  order  to  create  a  more  unified  set  of  detection  criteria,  and 
in  order  to  reach  fainter  sources.  X-ray  data  for  the  nine  smaller- 
field  clusters  with  X-ray  data  were  reprocessed  to  achieve  inter¬ 
nal  consistency  and  maximize  source  detection.  We  used  Chan¬ 
dra  pipeline  DS  8.4.5.  Detailed  X-ray  data  and  analysis  have 
already  been  published  for  the  very  deep  observation  of  the 
ONC  (Feigelson  et  al.  2005;  Getman  et  al.  2005),  and  for  two 
shallower  fields  north  and  south  of  the  ONC  (Ramirez  et  al. 
2004b);  similarly,  for  NGC  2264,  Flaccomio  et  al.  (2006)  and 
Ramirez  et  al.  (2004a)  have  published  deep  Chandra  data  (our 
region  here  is  primarily  covered  by  Flaccomio  et  al.  2006). 
Because  those  two  regions  include  X-ray  data  substantially 
deeper  than  those  for  the  rest  of  the  clusters,  we  have  not  re¬ 
processed  these  X-ray  data  here,  but  instead  taken  those  source 
lists  from  the  literature.  (A  primary  reason  we  reprocessed  the 
nine  smaller-field  clusters  is  to  reach  fainter  sources,  which 
the  deeper  integrations  in  Orion  and  NGC  2264  already  ac¬ 
complish.)  Additional  X-ray  data  obtained  contemporaneously 
with  Spilzer  monitoring  will  be  discussed  in  the  appropriate 
cluster-specific  papers. 

The  region  with  ACIS-I  data  is  usually  smaller  than  the  full 
region  observed  with  Spitzer  during  the  cryogenic  era,  but  often 
covers  the  region  monitored  for  YSOVAR;  see  Figures  2-13. 
Aside  from  AFGL  490  (which  has  no  X-ray  data),  the  region 
with  the  least  spatial  Chandra  fractional  coverage  is  Orion;  the 
Orion  region  with  IR  light  curves  is  by  far  the  largest  of  our 
sample.  The  other  10  clusters  have  complete  or  nearly  complete 
X-ray  coverage  of  the  region  with  two-band  IRAC  light  curves; 
X-ray  coverage  of  the  regions  with  light  curves  in  only  one- 
band  varies.  The  Chandra  sensitivity  varies  across  the  FOV, 
with  higher  sensitivity  in  the  center  of  the  pointing,  decreasing 
toward  the  edges. 

Source  detection  from  the  archival  Chandra  data  is  performed 
with  the  wavdetect  algorithm  in  CIAO  (Chandra  Interactive 
Analysis  of  Observations;  Fruscione  et  al.  2006),  versions  4.5 
and  4.6.  The  details  for  the  detection  process  used  here  are 
identical  to  those  used  in  Gunther  et  al.  (2012).  Our  chief 
goal  is  identification  of  all  X-ray  sources,  since  even  a  weak 
detection — provided  it  is  coincident  with  an  IRAC  source — has 
a  high  likelihood  of  being  a  legitimate  source,  though  this  does 
not  necessarily  prove  cluster  membership.  Table  6  lists  the 
number  of  X-ray  sources  detected  above  the  2 a  significance 
level  in  the  entire  ACIS-I  field  (a  single  ACIS-I  FOV  is 
~  1 6'  x  ~  1 6').  For  comparison.  Table  6  also  indicates  the  number 
of  detections  reported  in  the  literature,  often  for  exactly  the 
same  data.  We  recover  essentially  all  the  previously  reported 
sources,  and  add  a  significant  number  of  weaker  ones  due  to  the 
lower  threshold  employed  here.  Note  that  without  the  additional 
criterion  imposed  here  of  requiring  a  match  between  an  X-ray 
source  and  an  IRAC  source,  the  2 a  significance  level  used  here 
would  be  too  low  and  would  result  in  a  high  number  (about  one- 
third)  of  spurious  sources.  The  requirement  of  a  positional  match 
with  an  IRAC  source  allows  us  to  use  such  a  low  threshold. 
We  found  that  below  2a,  we  do  not  find  matches  between 
X-ray  and  IR  sources  in  excess  of  matches  expected  due  to 
random  chance. 

To  determine  if  a  given  X-ray  source  is  coincident  with  an 
IRAC  source,  we  matched  X-ray  sources  to  IRAC  source  lists 


generated  for  Spitzer  cryogenic-era  observations,  as  described 
in  Section  2.7  above.  Due  to  the  highly  spatially  dependent 
Chandra  point-spread  function,  three  matching  radii  were  used. 
For  sources  within  3'  of  the  aimpoint,  the  matching  radius  was 
1".  For  sources  more  than  6'  off-axis  from  the  aimpoint,  the 
matching  radius  was  2".  In  between,  the  matching  radius  was 
1"5.  We  note  that  IRAS  20050+2720  and  IC  1396A  are  both 
exceptions  since  the  fields  were  each  observed  multiple  times, 
at  different  roll  angles.  Since  IRAS  20050+2720  covered  a  wide 
range  of  rotation  angles,  in  that  case,  the  source  positions  are 
composites  of  different  point  spread  functions;  matches  were 
made  more  carefully  in  these  observations.  For  IC  1396A,  the 
range  of  roll  angles  was  smaller,  so  the  same  approach  was 
used  as  for  the  rest  of  the  clusters.  For  each  cluster,  in  Table  6, 
we  list  the  total  number  of  X-ray  sources  with  a  match  to  an 
IRAC  source  within  the  appropriate  radius.  We  note  here  for 
completeness  that  there  are  some  very  bright  X-ray  sources 
without  IR  counterparts;  for  our  purposes  within  YSOVAR,  we 
exclude  these  sources. 

If  an  X-ray  source  is  bright  enough  for  spectral  fitting,  one 
can  determine  four  key  characteristics;  X-ray  temperature,  gas 
absorption,  mean  flux  (Fx),  and  variability.  We  fit  a  single 
Astrophysical  Plasma  Emission  Code  (APEC)  thermal  emission 
model  (Smith  et  al.  2001)  for  each  source  detected  above  30 
net  counts,  using  C-statistics  and  leaving  the  absorbing  column 
density  (Ah),  the  temperature  (7),  and  the  volume  emission 
measure  as  free  parameters  to  determine  the  first  three  of 
the  above  characteristics.  The  metallicity  is  fixed  at  30%  of 
the  solar  abundances  (Anders  &  Grevesse  1989),  in  keeping 
with  typical  values  from  coronal  emission  from  late-type  stars. 
Details  of  the  extraction  and  fitting  processes  are  similar  to  those 
discussed  by  Winston  et  al.  (2010),  which  follows  the  procedure 
for  automated  processing  laid  out  for  the  ANCHORS  (AN 
archive  of  CHandra  Observations  of  Regions  of  Star  formation) 
pipeline  (Spitzbart  et  al.  2005).  The  key  point  is  that  all  fits 
are  done  assuming  that  the  source  is  a  star  with  a  thermal 
one-temperature  spectrum.  We  calculate  luminosities  (Lv)  for 
each  source  using  the  fluxes  (Fx)  from  0.3  keV  to  8.0  keV 
and  line-of-sight  absorptions,  assuming  the  distance  to  each 
cluster  given  in  Table  1.  Detailed  flux  errors  due  to  the  fit  were 
not  calculated  separately  for  each  source,  since  systematics  are 
likely  to  dominate;  instead,  global  flux  errors  were  determined 
using  the  CIAO  tool  dmextract.  This  process  calculates  a  simple 
error  estimate  based  on  photon  statistics  and  the  mean  value  of 
the  exposure  map  in  the  source  region.  These  errors  are  about 
4%  at  2000  counts,  about  35%  at  1 00  counts,  and  the  errors  reach 
100%  below  about  50  counts.  The  error  budget  is  dominated  by 
photon  counts  and  uncertainty  in  Ah-  There  is  no  evidence  that 
such  errors  are  markedly  biased  by  the  C-statistic. 

Since  errors  on  Lx  can  be  dominated  by  systematic  effects,  we 
have  reprocessed  all  sources  (in  those  nine  smaller-field  clusters 
with  Chandra  data),  even  those  with  previously  published 
fluxes,  to  ensure  uniform  spectral  fitting  methodology.  For  faint 
sources  with  fewer  than  30  counts,  no  fit  can  be  performed. 
In  this  case,  we  determine  a  median  photon  energy,  which, 
when  combined  with  the  count  rate,  leads  to  an  approximate 
flux  determination.  All  spectral  properties  presented  here  are 
effectively  time  averaged  over  all  observations. 

To  determine  the  fourth  of  the  characteristics  above  (vari¬ 
ability),  we  tested  the  light  curves  for  variability  using  the 
Gregory-Loredo  method  (GL-vary;  Gregory  &  Loredo  1992). 
This  method  uses  maximum-likelihood  statistics  and  evaluates 
a  large  number  of  possible  break  points  from  the  prediction  of 
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Table  5 

Summary  of  X-ray  Observations 


Cluster 

Aimpoint a  (J2000) 

Obsid 

Exp.  Time  (ks) 

Obs  Date 

Notes 

AFGL  490 

No  X-ray  observations 

NGC  1333 

03:29:05.60,  +31:19:19.00 

642 

43.2 

2000 Jul  12 

See  also  analysis  by  Getman  et  al.  (2002), 
Winston  et  al.  (2010) 

NGC  1333 

03:29:02.00,  +31:20:54.00 

6436 

39.5 

2006  Jul  5 

(As  above) 

NGC  1333 

03:29:02.00,  +31:20:54.00 

6437 

36.6 

2006  Jul  1 1 

(As  above) 

Orion 

Reanalysis  beyond  the  scope  of  this  paper; 
values  taken  from  Getman  et  al.  (2005)  and 
Ramirez  et  al.  (2004b) 

Mon  R2 

06:07:49.50,  -06:22:54.70 

1882 

98.1 

2000  Dec  2 

See  also  analysis  by  Kohno  et  al.  (2002), 
Nakajima  et  al.  (2003) 

GGD  12-15 

06  10  50.00,  -06  12  00.00 

12392 

67.3 

2011  Dec  23 

X-ray  catalog  not  yet  in  the  literature;  X-ray 
properties  to  be  discussed  in  detail  in  the 
corresponding  YSOVAR  cluster  paper 

NGC  2264 

Reanalysis  beyond  the  scope  of  this  paper; 
values  taken  from  Flaccomio  et  al.  (2006) 
and  Ramirez  et  al.  (2004a) 

L1688 

16  27  17.18,  -24  34  39.00 

635 

100.7 

2001  May  16 

See  also  analysis  by  Imanishi  et  al.  (2001) 

Serpens  Main 

18  29  50.00, +01  15  30.00 

4479 

88.5 

2005  Jun  28 

See  also  analysis  by  Winston  et  al.  (2007), 
Giardino  et  al.  (2007) 

Serpens  South 

18  30  03.00,  -02  01  58.20 

11013 

99.5 

2011 Jun  10 

X-ray  catalog  not  yet  in  the  literature;  X-ray 
properties  to  be  discussed  in  detail  in  the 
corresponding  YSOVAR  cluster  paper 

IRAS  20050+2720 

20  07  13.60, +27  28  48.80 

6438 

22.6 

2006  Dec  10 

See  also  analysis  by  Giinther  et  al.  (2012) 

IRAS  20050+2720 

20  07  13.60,  +27  28  48.80 

7254 

20.9 

2007  Jun  7 

(As  above) 

IRAS  20050+2720 

20  07  13.60,  +27  28  48.80 

8492 

.  50.5 

2007  Jan  29 

(As  above) 

IC1396A 

21  36  50.30, +57  30  24.00 

10990 

29.7 

2010 Jun 9 

See  also  analysis  by  in  Getman  et  al.  (2012) 

1C1396A 

21  36  50.30,  +57  30  24.00 

11807 

29.8 

2010  Mar  31 

(As  above) 

CephC 

23  05  51.00,  +62  30  55.00 

10934 

44 

2010  Sep  21 

X-ray  catalog  not  yet  in  the  literature;  X-ray 
properties  to  be  discussed  in  detail  in  the 
corresponding  YSOVAR  cluster  paper 

Note.  a  Formally,  “aimpoint"  is  the  location  in  chip  coordinates  on  the  Chandra  detector  where  the  target  lands.  This  is  different  from  the  focal  point,  which  is  the 
location  of  the  sharpest  (narrowest)  point-spread  function.  In  practice,  this  can  be  thought  of  as  the  target  of  the  observation. 


Table  6 

X-ray  Source  Detection  Characteristics 


Cluster 

No.  >2era 

No.  Cryo  Match*1 

Detections  in  Literature 

Min(log  Lxf 

Maxllog  Lxf 

Medflog  Lx)c 

AFGL  490 

(No  X-ray  data) 

NGC  1333 

192 

119 

109  (Getman  et  al.  2002),  193  (Winston  et  al.  2010) 

27.94 

30.74 

29.12 

Orion 

Using  lit.;  see  text 

27.58 

32.62 

29.72 

Mon  R2 

492 

167 

154  (5rr;  Kohno  el  al.  2002) 

29.58 

31.16 

30.16 

GGD  12-15 

229 

172 

29.66 

31.77 

30.33 

NGC  2264 

Using  lit.;  see  text 

28,45 

31.29 

29.95 

L 1 688 

315 

69 

1 1  (Imanishi  et  al.  200 1 ) 

27.48 

31.30 

29.51 

Serpens  Main 

204 

161 

85  (Giardino  et  al.  2007) 

29.13 

31.71 

29.97 

Serpens  South 

294 

82 

29.00 

31.27 

29.82 

IRAS  20050+2720 

348 

239 

239  (Gunther  et  al.  2012) 

29.28 

30.88 

30.12 

1C  1396A 

185 

129 

415d  (Getman  et  al.  2012) 

29.14 

32.26 

30.34 

Ceph  C 

200 

97 

29.44 

30.95 

30.17 

Notes. 

a  Number  X-ray  sources  detected  via  our  reprocessing  at  >2tr. 
b  Number  of  our  X-ray  sources  matched  to  sources  in  the  1RAC  cryogenic-era  catalog. 

c  For  the  sample  of  X-ray  detected  sources  with  YSOVAR  light  curves  (the  subset  of  the  standard  sample  for  statistics  detected  in  X-rays),  the  minimum,  maximum, 
and  median  log  Lx,  where  L,  is  in  erg  s“ 1 . 

d  The  detections  used  in  Getman  et  al.  (2012)  do  not  require  1R  matches,  and  they  go  well  below  our  typical  significance  cut  off. 


constancy.  It  assigns  an  index  to  each  light  curve — the  higher 
the  value  of  the  index,  the  greater  the  variability.  Index  values 
greater  than  7  indicate  >99%  variability  probability.  Values  of 
the  GL-vary  index  >9  usually  indicate  flares.  GL-vary  is  not 
reliable  below  about  30  raw  counts,  the  same  limit  we  used  for 
performing  spectral  fits.  The  value  of  this  index  will  be  provided 


where  relevant  on  a  source-by-source  basis  in  the  individual 
cluster  papers. 

There  are  five  broad  classes  of  X-ray  sources  in  the  field. 

1 .  First,  background  active  galactic  nuclei  (AGNs) — these  will 
be  numerous;  up  to  50  of  these  per  16'  x  16'  Chandra  field 


22 


The  Astronomical  Journal,  148:92  (46pp),  2014  November 


Rebull  et  al. 


are  expected,  depending  on  the  depth  of  exposure  (Getman 
et  al,  2006).  However,  they  tend  to  be  faint  in  both  X-rays 
and  mid-IR,  or  are  not  matched  in  the  IRAC  bands  at  all. 

2.  Second,  background  starburst  galaxies — while  much  less 
common  than  AGNs,  they  are  brighter  than  AGNs,  and  have 
colors  similar  to  Class  II  YSOs.  They  tend  to  be  fainter  in 
the  IR  than  typical  Class  II  YSOs,  and  may  be  tentatively 
identified  via  the  faintness  of  the  IR  counterpart. 

3.  Third,  compact  objects  such  as  white  dwarfs — these  tend 
to  be  very  faint  and  usually  undetected  in  the  IRAC  bands. 

4.  Fourth,  YSOs — those  with  disks  have  already  been  identi¬ 
fied  via  their  IR  excesses.  We  identify  the  probable  disk- 
free  objects  that  are  detected  in  X-rays  by  their  star-like 
IR  colors  and  magnitudes  consistent  with  membership.  In 
addition  to  matching  an  IRAC  source,  the  star-like  colors 
are  required  to  ensure  that  any  newly  revealed  sources  are 
probable  Class  III  objects  and  not  distant  starburst  galaxies. 

5.  Fifth  and  finally,  active  late-type  field  stars — both  fore¬ 
ground  and  background  stars  can  appear  with  X-ray  fluxes 
comparable  to  our  targets.  Based  on  a  study  of  IC  1396, 
Getman  et  al.  (2006)  estimate  that  there  could  be  <10  of 
these  per  16'  x  16'  Chandra  field.  Since  the  contaminants 
are  a  mixture  of  foreground  and  background  objects  of 
comparable  fluxes  to  our  targets  (with  less  foreground  and 
more  background  contamination  likely  for  closer  clusters), 
without  optical  spectra,  they  are  very  hard  to  discern  from 
Class  III  objects.  It  is  necessary,  then,  that  these  remain  in 
our  sample  of  candidate  members  selected  via  X-rays  and 
represent  a  source  of  contamination. 

We  note  that  Getman  et  al.  (2012)  estimate  (again  for 
IC  1396,  a  region  in  the  Galactic  plane)  a  22%  probability  that 
any  X-ray  source  with  any  IRAC  counterpart  is  not  a  member 
of  the  cluster.  That  rate  drops  to  about  10%  in  IC  1396  if  one 
eliminates  sources  without  2MASS  JHKS  detections.  Because 
the  contamination  rate  is  affected  by  absorption,  exposure  time, 
and  the  depth  to  which  one  extracts  sources,  it  may  be  different 
for  the  other  clusters. 

We  can  improve  our  inventory  of  YSOs  in  these  clusters  by 
identifying  objects  with  X-ray  detections,  IRAC  counterparts, 
and  SEDs  that  are  consistent  with  those  of  stars.  We  have 
adopted  this  YSO  selection  mechanism  as  the  other  main 
component  of  the  primary  sample  definitions;  see  Section  3. 
Note  that  we  define  SEDs  consistent  with  stars  to  be  those  with 
a  fitted  SED  Class  III  (see  Appendix  B),  but  that  there  is  room 
to  create  an  augmented  membership  list  (Section  3.3)  to  include 
objects  that  have  an  X-ray  detection,  an  IRAC  counterpart,  were 
not  identified  as  a  YSO  from  the  IR  alone,  but  that  have  an  SED 
consistent  with  a  YSO.  With  regard  to  foreground  or  background 
stellar  contamination,  because  any  appropriate  brightness  cutoff 
is  a  function  of  cluster  distance  (and  Ay),  we  defer  any  detailed 
exclusion  of  likely  foreground  or  background  stars  from  the 
member  sample  to  the  individual  cluster  papers. 

3.  SAMPLE  DEFINITIONS 

For  most  of  our  clusters,  there  are  no  well-established  mem¬ 
bership  lists;  the  exceptions  are  Orion  and  NGC  2264.  More¬ 
over,  the  membership  lists  in  the  literature  use  a  wide  variety 
of  wavelengths  and  survey  depths  to  identify  members,  and  the 
spectroscopic  follow-up  of  candidate  members  is  uneven.  If  we 
decided  to  depend  on  the  robust  member  identifications  only 
from  the  literature,  our  sample  would  be  greatly  reduced  for 


most  clusters  and  highly  biased.  Certainly,  our  variability  sur¬ 
vey  will  identify  new  candidate  members  based  on  the  light 
curve  properties.  To  attempt  to  make  fair  comparisons  between 
clusters,  we  need  to  define  a  set  of  (candidate)  members  in  the 
same  or  at  least  consistent  ways  between  clusters. 

As  discussed  above,  even  within  our  survey,  different  clusters 
may  have  different  amounts  of  monitoring  data  beyond  the 
“fast  cadence”  data.  Thus,  for  making  comparisons  between 
clusters,  we  need  to  define  a  standard  set  of  data  that  are  used 
for  calculating  statistics  and  identifying  variables. 

Therefore,  we  define  a  “standard  YSOVAR  sample,”  which  is 
the  sample  that  is  (primarily)  discussed  in  this  paper  and  which 
forms  the  common  core  of  the  papers  planned  for  each  cluster. 
Each  cluster  may  have  an  additional  “augmented  sample”  as 
well,  to  take  into  account  additional  member  identifications  from 
the  literature  or  our  own  data  where  possible  and  necessary.  We 
now  discuss  the  definitions  of  these  samples. 

3.1.  Standard  Set  of  Members 

Each  cluster  has  an  IR-selected  sample  of  member  candidates 
defined  by  the  Gutermuth  et  al.  (2008b,  2010)  and  G09  selec¬ 
tion  algorithm,  run  on  the  cryo-era  catalogs  created  anew  as  per 
the  methodology  above  (Section  2.7).  A  detailed  description  of 
the  YSO  candidate  color  selection  algorithm  can  be  found  in 
Appendix  A  of  G09.  A  sample  of  YSOs  selected  this  way  is 
thought  to  be  a  statistically  well-defined  sample  (see  statisti¬ 
cal  discussions  in  Gutermuth  et  al.  2008b  and  G09),  comprised 
almost  entirely  of  members,  though  some  contamination  from 
background  galaxies  or  asymptotic  giant  branch  stars  is  always 
possible.  Very  few  of  these  IR-selected  members  have  spec¬ 
troscopic  follow-up  (or  spectra  in  the  literature  pre-dating  the 
Spitzer  observations),  and,  as  such,  many  should  technically  be 
thought  of  as  YSO  candidates,  though  we  include  them  all  in 
the  set  of  members. 

All  clusters  except  AFGL  490  have  X-ray  data,  and  thus  also 
an  X-ray-selected  sample  of  YSOs.  (As  with  the  IR-selected 
sample,  few  of  these  have  spectra,  so  technically  they  are  YSO 
candidates.)  Since  we  should  have  identified  most  of  the  YSOs 
with  disks  (at  least  disks  detectable  at  3.6-24 /rm)  in  the  IR 
selection  process,  we  use  the  X-ray  data  to  identify  additional 
young  stars  without  disks  (e.g.,  Class  Ills).  Thus,  we  add  to  the 
set  of  IR-selected  members  the  X-ray-selected  sample  defined 
by  the  algorithm  described  above  in  Section  2. 10,  which  can  be 
summarized  as  objects  having  an  X-ray  detection  above  a  2cr 
significance  threshold,  having  a  match  to  an  IRAC  object  in  the 
cryo-era  Spitzer  catalog,  and  having  an  SED  shape  consistent 
with  it  being  a  disk-free  YSO  or  a  star  (e.g.,  Class  III;  see 
Appendix  B.2).  This  sample  is  specifically  designed  to  find 
and  add  to  our  set  of  members  those  members  without  disks. 
However,  it  should  be  noted  that  (1)  the  Galactic  contamination 
rate  is  likely  to  be  higher  in  this  X-ray-selected  sample  than 
in  the  IR-selected  sample;  and  (2)  specifically  because  of 
the  contamination  rates,  members  with  disks  are  identified  as 
members  from  the  IR  excess,  not  the  X-ray  flux,  though,  of 
course,  members  with  disks  can  also  have  measured  X-ray  fluxes 
in  our  database. 

We  have  thus  defined  our  “standard  set  of  members”  to 
be  the  union  of  all  IR-selected  members  with  disks  and 
X-ray-selected  members  without  disks.  There  are  provisions 
for  adding  additional  objects;  see  Section  3.3  below.  Note  that 
this  definition  can  be  applied  independently  of  whether  or  not 
there  is  a  light  curve,  but  of  course  in  the  context  of  this  dis¬ 
cussion  of  YSOVAR  data  we  require  a  light  curve.  Note  also 
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that  the  IR  selection  requires  four  bands  of  IRAC,  and  thus  both 
very  faint  and  very  bright  previously  identified  members  may 
be  omitted  from  the  standard  set;  objects  such  as  these  known 
to  be  members  via  some  other  approach  in  the  literature  may  be 
added  in  the  augmented  sample  (Section  3.3).  Finally,  note  that 
because  the  data  that  go  into  our  selection  of  cluster  members 
are  of  various  depths,  and  because  the  clusters  are  at  a  variety 
of  distances,  the  effective  mass  limit  reached  by  each  set  of 
standard  members  varies  from  cluster  to  cluster. 

3.2.  Standard  Set  for  Statistics 

All  of  the  original  YSOVAR  light  curves  were  obtained  with 
very  similar  HDR  mapping  observations  (see  Section  2.2)  in 
a  fast  cadence  (see  Section  2.4),  though  the  length  varies,  and 
some  clusters  have  additional  slow  cadence  observations  and/or 
staring  observations.  The  time  sampling  and  total  length  of  light 
curves  obtained  within  a  single  cluster  field  can  also  vary  as 
the  FOV  changes  with  time  (Section  2.3).  We  have  attempted 
to  remove  all  instrumental  effects  from  the  input  data,  and  only 
retained  photometry  where  there  were  valid  measurements  on 
at  least  three  BCD  frames  (Section  2.5). 

We  now  define  the  “standard  set  of  data  for  statistics”  as 
follows.  Since  the  fast  cadence  is  the  most  common  (and  most 
similar)  among  the  YSOVAR  clusters,  we  used  only  these 
mapping  fast  cadence  data,  for  those  light  curves  that  have 
at  least  five  viable  epochs  (with  each  epoch  obtained  from  at 
least  three  BCDs  per  epoch  that  are  not  obviously  compromised 
by  instrumental  effects  or  cosmic  rays),  to  calculate  statistical 
quantities  such  as  mean,  median,  etc.,  as  well  as  Stetson  index 
and  x2  (discussed  below;  see  Section  5).  We  have  defined 
statistical  values  calculated  on  the  fast  cadence  data  as  the 
“standard  set  of  statistical  values”  for  each  cluster  and  employ 
them  to  identify  variables  and  compare  values  across  clusters. 
Finally,  for  stars  fainter  than  [3.6]— [4.5]  ~  16,  noise  tends 
to  dominate  the  light  curves  (Sections  4.2  and  5).  We  have 
retained  these  faint  sources,  but  objects  this  faint  are  considered 
individually  where  relevant. 

Note  that  the  standard  set  for  statistics  is  thus  defined  as  all 
light  curves  with  at  least  five  points,  just  in  the  fast  cadence. 
This  is  independent  of  whether  or  not  the  target  is  identified  as 
a  member. 

Elsewhere  in  this  paper,  we  refer  to  “all  objects  with  a  light 
curve” — this  means  anything  with  a  light  curve  in  the  standard 
set  for  statistics.  (Essentially  no  YSOVAR-classic  sources  have 
only  points  outside  of  the  fast  cadence.) 

When  available,  additional  epochs  of  data  can  be  included 
for  additional  calculations  on  a  per-cluster  basis  in  the  corre¬ 
sponding  papers,  and  will  clearly  be  indicated  as  such  where 
relevant. 

3.3.  Augmented  Sample  of  Members 

We  identified  members  above  using  IR  and  X-ray  data,  which 
implicitly  relies  on  the  shape  of  the  SED  between  2  and  24  pm. 
That  sample  is  still  the  best  set  of  members  that  we  will  use 
to  compare  across  clusters,  because  that  set  of  members  is 
defined  as  similarly  as  possible  across  all  clusters.  However, 
additional  young  objects  may  be  identified  in  the  literature,  and 
additional  members  may  be  suggested  based  on  our  own  data. 
The  “augmented  sample  of  members”  is  where  these  additional 
likely  members  can  be  included. 

Some  clusters  (e.g.,  NGC  1333,  LI 688)  have  considerable 
literature  discussion  of  members,  and  others  (e.g.,  AFGL  490, 


Mon  R2)  have  far  less.  We  therefore  cannot  rely  exclusively  on 
the  literature  to  select  members,  but  neither  should  we  ignore 
members  identified  in  the  literature  and  not  selected  above. 
Thus,  each  cluster  paper  may  include  in  the  augmented  sample 
the  literature-identified  sources  that  are  not  already  found  using 
our  IR  or  X-ray  methods  above. 

We  can  use  our  own  data  to  identify  new  cluster  members. 
While  only  1  %-2%  of  the  field  population  may  be  variable,  YSO 
variability  is  the  rule,  not  the  exception.  It  is  therefore  possible  to 
identify  new  cluster  members  from  light  curve  properties  alone; 
one  could  identify  all  variables  as  new  members,  or  one  could 
take  only  those  with  certain  properties  such  as  amplitude  above 
a  threshold.  We  could  identify  cluster  members  from  either  the 
standard  set  for  statistics  (the  fast  cadence),  or,  for  those  clusters 
with  longer  cadences,  from  those  additional  data. 

The  set  of  statistically  selected  variables  (see  Section  5) 
are  those  with  Stetson  index  greater  than  0.9,  and/or  with  /2 
greater  than  5,  and/or  with  a  significant  period,  calculated  over 
only  the  standard  set  for  statistics  (the  fast  cadence).  Often, 
these  variables  should  also  be  identified  as  cluster  members, 
but  these  individual  objects  will  be  discussed  on  a  per-cluster 
basis  (because,  for  example,  they  can  be  background  eclipsing 
binaries;  Morales-Calderon  et  al.  2012).  Variables  identified 
using  data  beyond  the  standard  set  for  statistics  (data  beyond 
the  fast  cadence  data)  will  also  be  included  on  a  per-cluster 
basis,  and  will  be  discussed  in  the  cluster  papers.  Those 
newly  identified  cluster  members  may  also  be  included  in  the 
augmented  sample  of  members. 

While  variability-identified  objects  will  make  an  important 
contribution  to  our  understanding  of  the  complete  membership 
of  each  cluster,  they  should  not  be  used  in  calculations  of,  e.g., 
variability  fractions;  that  sample  should  be  selected  on  the  basis 
of  a  parameter  distinct  from  variability,  such  as  disk  excess  or 
X-ray  emission.  This  is  why  the  new  candidate  members  we 
identify  from  our  data  are  included  in  the  augmented  sample  of 
members  and  not  in  the  standard  set  of  members. 

This  augmented  set  of  members  is  only  used  (where  it  is 
clearly  identified)  in  the  individual  cluster  papers,  not  in  the 
remainder  of  this  paper. 

4.  ENSEMBLE  ANALYSIS 

In  this  section,  we  present  an  analysis  of  the  entire  set  of  data 
we  used  for  our  clusters,  independent  of  variability,  which  is 
discussed  in  Section  5. 

4.1.  Cluster  Parameterization 

In  order  to  compare  results  among  clusters,  it  is  useful  to 
be  able  to  place  clusters  in  some  sort  of  relative  order  that 
could,  in  the  most  useful  (though  hypothetical)  case,  be  tied  to 
age.  There  are  various  ways  of  parameterizing  the  evolutionary 
state  of  these  clusters,  and  we  considered  several,  all  aimed  at 
capturing  the  relative  numbers  of  sources  in  various  SED  class 
bins  (see  Appendix  B  for  a  brief  definition  of  SED  classes  and 
our  placement  of  sources  therein).  Such  ratios  in  some  sense 
capture  the  relative  “degree  of  embeddedness”  of  sources  in 
these  regions,  perhaps  with  an  ultimate  (though  undefined  here) 
link  to  cluster  ages.  Formally,  we  are  binning  by  SED  slope, 
so  the  parameterizations  are,  strictly  speaking,  relative  fractions 
of  sources  with  a  given  SED  slope,  indicative  of  the  amount  of 
circumstellar  material,  e.g.,  how  self-embedded  a  given  source 
may  be.  We  interpret  clusters  with  more  sources  with  large 
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SED  slopes  to,  on  average,  contain  more  sources  that  are  more 
embedded. 

G09  chose  to  use  the  ratio  of  Class  II  to  Class  I  sources.  These 
ratios  were  all  obtained  internally  consistently,  e.g.,  sources 
selected  and  categorized  according  to  the  same  series  of  IR 
color  cuts  and  data  reduction.  Because  G09  was  working  with 
Spitz.er  data  over  a  relatively  large  region  for  each  cluster,  this 
Class  II  to  Class  I  ratio  could  be  calculated  for  the  whole  region 
and  for  subsections  of  the  region.  The  Class  II  to  Class  I  ratios 
as  calculated  for  the  cores  of  these  clusters  (e.g..  Table  6  in  G09) 
appear  in  our  Table  1  and  in  the  discussion  in  Section  2. 1 .2. 

In  YSOVAR,  we  have  clusters  not  included  in  G09,  and 
moreover,  even  for  the  clusters  included  in  G09,  we  generally 
have  light  curves  for  only  a  small  subset  of  the  sources  G09 
considered  because  we  observed  a  smaller  region.  To  calculate 
a  Class  II  to  Class  I  ratio  for  the  portion  of  each  cluster 
sampled  by  the  YSOVAR  monitoring,  we  performed  the  same 
calculation  for  objects  in  the  standard  set  for  statistics  (e.g., 
having  YSOVAR  light  curves)  by  reducing  the  cryogenic  data 
in  the  same  way  and  performing  the  same  series  of  G09  color 
cuts  and  classification;  see  Section  2.7.  Then,  we  recalculated 
the  ratio  of  Class  II  to  Class  I  sources  specific  to  the  YSOVAR 
data  using  the  classes  assigned  via  the  G09  algorithm  only  for 
those  sources  with  light  curves.  Those  Class  II  to  Class  I  ratios 
also  appear  in  Table  I  and  in  the  discussion  in  Section  2.1.2. 
The  values  are  provided  in  the  present  work  in  part  as  a  link 
to  the  G09  analysis;  note  that  they  are  calculated  in  the  same 
way  as  G09,  e.g.,  with  the  IR-selected  sources  alone,  not  on  the 
standard  set  of  members  per  se. 

The  G09  parameterization  is  based  only  on  IR-selected 
sources.  For  most  of  our  clusters,  we  have  X-ray  data  as  well 
(see  Section  2.10),  so  we  at  least  have  some  information  on  the 
Class  III  population.  It  is,  however,  true  that  we  do  not  always 
have  complete  Chandra  coverage  of  our  fields  (further  discus¬ 
sion  of  coverage  appears  in  Section  2.10  and  Figures  1-13), 
and  even  for  clusters  where  we  have  Chandra  coverage,  the 
Chandra  sensitivity  is  a  strong  function  of  location  on  the 
array.  Nonetheless,  we  would  like  to  include  the  information 
we  have,  and  simply  using  the  Class  II  to  Class  I  ratio  does  not 
incorporate  information  about  the  Class  III  population. 

We  explored  several  alternative  parameterizations  of  the 
relative  fractions  of  embedded  sources,  all  of  which  involved 
various  ratios  of  classes  (or  groups  of  classes)  to  the  total  or 
other  classes  (or  groups  of  classes).  Histograms  of  the  relative 
fractions  of  the  SED  classes  for  the  standard  set  of  members 
(Section  3. 1 )  for  objects  with  light  curves  in  the  standard  set  for 
statistics  (Section  3.2)  in  each  cluster  appear  in  Figure  1 7.  AFGL 
490  can  be  seen  to  be  deficient  in  a  complete  sample  of  Class  III 
objects  because  it  has  no  X-ray  data;  the  Class  III  objects 
identified  here  have  small  IR  excesses  (or  sufficient  reddening 
at  2 /urn)  and  thus  SED  slopes  consistent  with  Class  III.  By 
many  metrics,  Serpens  South  has  the  highest  fraction  of  sources 
with  positive  SED  slopes,  which  we  take  to  be  most  embedded 
sources.  Similarly,  both  Orion  and  IC  1396A  have  the  highest 
fraction  of  sources  with  negative  SED  slopes,  which  we  take  to 
be  less  embedded  sources. 

For  further  analysis  here,  we  have  settled  on  the  fraction  of 
Class  I  sources  to  the  total  number  of  members,  for  objects  with 
light  curves  (objects  in  the  standard  set  for  statistics).  These 
ratios  are  also  included  in  Table  1 .  They  include  (as  part  of  the 
total  number  of  members)  objects  selected  via  X-rays.  However, 
there  are  still  the  fundamental,  systematic  uncertainties  inherent 
in  the  classification  approach,  in  the  selection  of  members 


class  class 

Figure  17.  Histograms  of  the  relative  fractions  of  fitted  SED  classes  (classes 
derived  from  SED  fits  as  discussed  in  Appendix  B.2)  for  the  standard  set  of 
members  with  light  curves.  Clusters  appear  in  R.A.  order.  By  nearly  any  ratio  of 
classes  used  as  a  parameterization,  IC  1 396A  and  Orion  have  the  highest  fraction 
of  sources  with  negative  SED  slopes  (taken  to  be  less  embedded  sources),  and 
Serpens  South  has  the  highest  fraction  of  sources  with  positive  SED  slopes 
(taken  to  be  most  embedded  sources).  AFGL  490  is  noticeably  incomplete 
in  the  Class  III  bin;  there  are  no  X-ray  observations  available  for  it,  and  the 
Class  III  objects  in  that  bin  have  small  IR  excesses  and  SED  slopes  consistent 
with  Class  III  (see  Appendix  B  for  a  description  of  classes  and  class  selection). 

without  spectroscopic  follow-up,  in  the  completeness  of  the 
surveys  involved  (in  both  area  and  depth),  and  the  requirement 
that  there  be  a  light  curve  (e.g.,  bright  enough  in  the  IRAC 
channels)  as  well  as  in  the  relative  paucity  of  Class  I  objects 
overall.  Additionally,  for  AFGL  490,  there  are  no  X-rays  to 
be  used,  so  the  ratio  is  calculated  with  solely  the  IR-identified 
members.  However,  the  Class  III  objects  appear  only  in  the 
denominator,  combined  with  all  the  other  classes.  We  also  note 
that  our  inventory  of  Class  I  sources  must  be  incomplete,  since 
we  lack  the  long  wavelength  coverage  that  would  be  needed 
to  find  the  most  embedded  sources  (e.g.,  Stutz  et  al.  2013),  but 
such  extremely  embedded  sources  will  also  not  have  a  YSOVAR 
light  curve. 

This  parameterization  using  the  Class  I/total  ratio  should  be 
related  to  the  G09  Class  II/Class  I  ratio  determined  for  the 
objects  having  light  curves.  Figure  18  plots  these  parameteri¬ 
zations  against  each  other.  The  error  bars  are  derived  assuming 
uncorrelated  Poisson  statistics,  because  it  is  difficult  to  quantify 
the  additional  systematic  errors  described  above.  The  best-fit 
slope  of  the  line  fit  to  this  relation  (taking  into  account  errors  in 
both  directions  on  each  point)  is  —0.024  ±  0.005.  The  correla¬ 
tion  coefficient,  Pearson’s  r,  calculated  for  these  two  parameters 
is  —0.87,  and  a  probability  that  the  parameterizations  are  not 
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Figure  18.  Comparison  of  the  Class  I  I/I  parameterization  from  G09  and  the  Class  I/total  parameterization  used  here.  All  values  are  calculated  only  for  objects  with 
light  curves  (the  standard  set  for  statistics;  Section  3.2);  the  G09  values  use  only  the  IR-selected  members  from  the  G09  approach,  and  the  Class  I/total  values  use 
the  standard  set  of  members  (Section  3.1).  The  error  bars  are  derived  assuming  uncorrelated  Poisson  statistics.  Cluster  labels  appear  at  the  Class  II/Class  I  location 
corresponding  to  the  cluster.  The  gray  line  is  the  best-fit  line,  using  errors  in  both  directions  on  each  point,  and  has  a  slope  of  —0.024  =k  0.005.  Pearson’s  correlation 
coefficient  (r),  calculated  for  these  points  is  —0.86;  the  calculated  probability  that  the  parameterizations  are  not  correlated  is  only  0.04%.  We  take  these  values  to  be 
correlated.  The  additional  gray  points  are  subregions  of  Orion,  used  here  to  show  the  scatter  inherent  in  the  large  Orion  map  (see  the  text).  Notional  “young”  and  “old” 
annotations  on  the  axes  describe  approximate  relative  ages  that  may  quantitatively  correspond  to  large  or  small  values  of  these  parameterizations. 


correlated  of  only  0.04%.  We  assume  based  on  the  statistics  and 
our  underlying  physical  intuition  that  these  values  are  indeed 
correlated. 

Several  individual  points  in  Figure  1 8  merit  additional  dis¬ 
cussion.  Despite  having  no  X-ray  data,  AFGL  490  is  consistent 
with  the  trend  shown  in  Figure  18.  (If  one  assumes  that  there 
might  be  about  as  many  Class  Ills  as  Class  IIs  in  this  region, 
then  the  point  could  move  down  to  about  0.1-0.15,  which  would 
still  be  broadly  consistent  with  the  trend.)  The  NGC  2264  point 
is  below  the  trend.  The  X-ray  data  obtained  for  NGC  2264  is 
deeper  than  the  X-ray  data  for  the  other  clusters,  except  for  the 
central  Orion  region.  This  results  in  more  of  the  fainter  sources 
being  included  in  the  total  number  of  YSOs,  and  pushes  the 
Class  I/total  ratio  toward  lower  numbers,  as  seen.  For  both  pa¬ 
rameterizations,  Serpens  South  is  selected  as  the  cluster  with 
the  most  embedded  sources,  as  expected  from  Figure  17.  It  is 
well  above  the  fitted  line  in  Figure  18;  perhaps  an  exponential 
decay  rather  than  a  simple  line  would  be  a  better  fit  to  use,  but, 
in  the  absence  of  additional  very  embedded  clusters  to  constrain 
the  most  embedded  end  of  the  distribution  (or,  indeed,  spectro¬ 
scopic  vetting  of  the  members  and  the  reduction  of  other  such 
uncertainties),  a  line  is  the  simplest  fit  to  use.  IC  1396A.  based 
on  Figure  17,  should  be  one  of  the  clusters  with  the  fewest  em¬ 
bedded  sources.  It  is  identified  as  the  least  embedded  using  the 
Class  I/total  ratio;  there  is  a  large  uncertainty  on  the  Class  II/ 
Class  I  ratio,  and  it  is  consistent  with  being  the  least  embedded 
within  lcr.  The  Class  II/Class  I  ratio  formally  identifies  Orion 


as  the  least  embedded.  However,  aside  from  AFGL  490  where 
there  are  no  X-ray  data,  Orion  is  the  cluster  with  the  poorest 
match  between  the  YSOVAR-monitored  region  and  the  existing 
X-ray  data  coverage.  Orion  has  a  very  deep  X-ray  pointing,  but 
only  in  the  central  ONC  region  (Feigelson  et  al.  2005;  Getman 
et  al.  2005).  There  are  two  shallower  pointings  that  contribute  X- 
ray  data  (see  Section  2.10  and  Figure  3),  but  the  region  of  sky  in 
Orion  for  which  we  have  IR  light  curves  has  the  least  fractional 
coverage  in  X-rays  of  all  our  clusters  (aside  from  AFGL  490). 
To  investigate  the  degree  to  which  the  uneven  X-ray  coverage 
affects  the  placement  of  Orion  in  this  diagram,  Figure  1 8  also  in¬ 
cludes  points  for  Orion  when  broken  into  “North”  (declination  > 
—05:05:25),  “South”  (declination  <—05:33:15°),  and  “ONC” 
(between  those  two  limits)  fields.  There  is  scatter,  clearly,  in 
these  points,  but  when  those  points  are  used  instead  of  the  sin¬ 
gle  Orion  point,  the  fit  is  functionally  indistinguishable  from 
that  using  the  single  Orion  point. 

We  use  the  Class  I/total  parameterization  in  subsequent 
discussions  in  this  paper,  most  notably  Section  6. 

4.2.  Brightness  Distribution  at  J,  [3.6],  and  [4.5] 

Figures  19  and  20  show  the  distribution  of  J,  [3.6],  and 
[4.5],  in  units  of  magnitude  and  the  percent  of  the  sample, 
for  the  standard  set  for  statistics  (all  objects  with  at  least 
five  points  in  the  YSOVAR  fast-cadence  light  curve,  including 
cluster  members  and  field  objects).  The  total  number  of  objects 
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Figure  19.  Distribution  of  J  magnitudes  (left)  and  IRAC  magnitudes  (right;  [3.6]  in  the  solid  line  and  [4.5J  in  the  dotted  line),  for  the  standard  set  for  statistics  (all 
objects  with  light  curves),  in  units  of  fraction  of  sample  (in  %)  for  each  cluster.  This  figure  shows  AFGL490,  NGC  1333,  Orion,  Mon  R2,  GGD  12-15,  and  NGC  2264. 
See  the  text  for  discussion.  As  a  result  of  these  plots,  we  are  cautious  about  objects  fainter  than  [3.6J— [4.5]  ~  16,  both  because  they  have  low  signal-to-noise  in  our 
monitoring  data  and  because  they  are  likely  dominated  by  non-members. 


portrayed  in  the  figures  ranges  from  ~  1 00  (J:  L 1 688  or  Serpens 
South)  to  ~7200  ([3.6]  and  [4.5]:  Orion).  The  J  data  are 
generally  from  2MASS;  notably,  NGC  1333  is  deeper  than 
the  other  clusters  in  J  because  additional  data  from  the  6  x 
2MASS  survey  have  been  included  (Section  2.9).  Generally, 
more  objects  are  available  at  [3.6]  or  [4.5]  than  at  J ,  largely  due 
to  the  greater  effective  depth  of  Spitzer.  In  several  cases  (most 
notably  LI 688,  Serpens  South,  and  NGC  2264),  a  relatively 
small  fraction  of  the  objects  with  Spitzer  light  curves  have  J 
counterparts,  since  these  clusters  are,  on  average,  generally  more 
embedded  than  the  others.  For  most  of  the  clusters,  for  most  of 
the  objects,  JHKS  data  are  not  available  for  objects  with  [3.6]  or 
[4.5]  fainter  than  about  15th  mag. 

The  Orion  maps  extend  out  beyond  the  edges  of  the  cluster, 
and  include  a  higher  proportion  of  field  stars  and  other  con¬ 
taminants  than  do  the  other  smaller-field  clusters.  This  can  be 
seen  in  the  structure  in  the  Orion  [3.6]  and  [4.5]  histogram, 
which  is  double-peaked;  the  brighter  peak  is  likely  dominated 
by  the  cluster  members,  and  the  fainter  peak  is  likely  domi¬ 
nated  by  contaminants.  For  similar  reasons,  if  the  UKIDSS  data 


(Section  2.9)  are  included,  the  Serpens  Main  J  histogram  ex¬ 
tends  to  J  ~  20  and  is  also  double-peaked,  but  the  LI 688  J 
histogram  is  not  so  obviously  double-peaked,  likely  due  to  the 
higher  obscuration  levels  of  the  background  population. 

Mon  R2  seems  to  be  different  from  the  other  clusters  in  that 
the  fainter  end  of  the  [3.6]  and  [4.5]  histograms  are  considerably 
flatter.  This  is  most  likely  a  symptom  of  the  difficulty  of 
obtaining  Spitzer  light  curves  for  faint  sources  in  the  presence 
of  high  and  spatially  variable  background;  there  are  some  very 
IR-bright  sources  in  the  IRAC  FOV  (see  Figure  4),  and  the 
scattered  light  is  substantial,  coupled  with  intrinsically  bright 
outflows  and  PAH  features.  No  UKIDSS  J  mag  in  this  region 
were  in  the  public  archive  at  the  time  we  checked  (in  2013 
September). 

As  can  be  seen  from  the  turnover  at  about  1 6th  mag  in  the 
[3.6]  and  [4.5]  histograms  in  Figures  19  and  20,  our  data  do  not 
extend  much  fainter  than  [3.6]-[4.5]  ~  16  mag.  By  inspection, 
all  of  these  faint  objects  appear  to  be  legitimate  point  sources  on 
the  images.  For  stars  fainter  than  this,  noise  tends  to  dominate 
the  light  curves  (see  Section  5).  It  is  hard  to  completely  reject 
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Figure  20.  Same  as  for  Figure  19,  but  for  L1688,  Serpens  Main,  Serpens  South,  IRAS  20050+2720,  IC  1396A,  and  Ceph  C. 


these  fainter  sources — we  could  exclude  those  whose  cryogenic- 
era  [3.6]  >  16,  though  the  YSOVAR  epochs  could  vary  above 
and  below  that  boundary;  or  we  could  discard  those  where  the 
mean  during  the  YSOVAR  campaign  is  >  16,  but  there  are  some 
objects  for  which  the  mean  [3.6]  <16  but  the  mean  [4.5]  >16. 
As  noted  in  Section  3,  we  have  retained  these  faint  sources,  but, 
in  general,  the  large  uncertainties  associated  with  their  Spitzer 
photometry  preclude  strong  statements  about  their  variability; 
objects  this  faint  are  considered  individually  where  relevant. 

Similar  histograms  for  only  the  standard  set  of  members 
(identified  through  IR  excess  and  X-ray  emission;  see  Section  3) 
are  considerably  less  populated,  and  brighter;  the  peaks  are 
around  [3 .6]— [4.5]  ~  12  mag.  This  is  consistent  with  the  location 
of  the  brighter  peak  in  Orion,  and  several  of  the  tails  of  the 
distributions  seen  in  Figures  19  and  20. 

4.3.  X-Ray  Brightness  Distribution 

Figure  21  contains  histograms  of  the  log  of  the  X-ray  lumi¬ 
nosities  (log  Lx,  where  Lx  is  in  erg  s-i)  for  those  objects  with 
light  curves  (standard  set  for  statistics)  and  bright  enough  in 
flux  (Fx)  to  have  a  calculated  Lx  (Section  2.10).  The  Orion 
and  NGC  2264  histograms  reach  fainter  values  in  Fx  than 
those  of  the  other  clusters  because  those  integrations  were 


considerably  deeper.  Moreover,  essentially  the  entire  NGC  2264 
field  considered  here  has  X-ray  data,  whereas  the  fractional 
X-ray  coverage  of  the  Orion  field  is  relatively  low  com¬ 
pared  to  NGC  2264  or  the  other  clusters  here  (aside  from 
AFGL  490,  where  there  is  no  X-ray  data).  Because  this  fig¬ 
ure  has  incorporated  distance  (distances  are  listed  in  Table  1 ) 
to  the  clusters  in  the  calculation  of  Lx,  the  two  closest  clus¬ 
ters  (NGC  1333  and  LI 688)  have  histograms  reaching  the 
faintest  Lx. 

5.  IDENTIFYING  VARIABLES 

There  are  many  ways  discussed  in  the  literature  of  identifying 
variables  in  time  series  data.  We  tested  several  methods,  and 
settled  on  three  primary  ones,  which  we  now  discuss  in  separate 
subsections.  Recall  that  we  are  calculating  statistics  for  the 
standard  statistical  sample,  e.g.,  on  only  the  fast  cadence  data, 
for  those  objects  with  at  least  five  viable  data  points  in  the  light 
curve  (Section  3). 

5.7.  Stetson  Index 

The  first  way  we  identify  variables  is  the  Stetson  index 
(Stetson  1996),  which  quantifies  correlation  of  variability  in 
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log(Lx)  log(Lx) 

Figure  21.  Distribution  of  log(Z-r)  values  (in  erg  s_l )  for  the  objects  in  the  standard  set  for  statistics  (with  Lx  detections)  in  our  clusters.  Note  that  there  are  no  data 
available  for  AFGL  490,  and  we  use  literature  data  for  NGC  2264  and  Orion;  see  the  text.  L1688  and  NGC  1333  are  the  closest  clusters,  and  so  relatively  faint  L , 
measurements  are  obtained  with  even  relatively  shallow  observations.  NGC  2264  and  Orion  have  the  deepest  integrations  and  therefore  also  include  relatively  faint  Lx. 


two  (or  more)  bands.  The  Stetson  variability  index  is  computed 
for  each  object  as: 


E,=i  8i  x  sgn(Pj)  x  y/\Pj\ 

Z-»i'= 1  Si 


(i) 


where  N  is  the  number  of  pairs  of  observations  for  a  star  taken  at 
the  same  time,33/3,  =  is  the  product  of  the  normalized 

residuals  of  two  observations,  and  g,  is  the  weight  assigned  to 
each  normalized  residual.  In  our  case,  the  weights  are  all  equal 
to  one.  The  normalized  residual  (5)  for  a  given  band  is  computed 
as: 

= 

V  N  —  1  ct, 

where  N  is  the  number  of  measurements  used  to  determine  the 
mean  magnitude  and  07  is  the  photometric  uncertainty.  Objects 


33  The  II  and  12  maps  are  taken  of  any  given  source  (providing  it  falls  within 
the  region  with  two-band  coverage)  typically  within  ^12  minutes  of  each 
other.  For  these  observations  taken  on  a  few-epochs-per-day  cadence,  we  are 
not  sensitive  to  timescales  of  minutes,  and  the  data  points  are  effectively 
simultaneous. 


with  larger  values  of  the  Stetson  index  are  typically  taken  to 
be  variable.  Since  errors  are  included  in  the  calculation,  light 
curves  that  are  just  noisy  are  not  identified  as  variable.  Objects 
with  variability  in  different  bands  that  is  not  correlated  will 
not  be  identified  via  this  method;  physically,  we  expect  most 
YSOs  to  have  similar  variations  in  the  two  IRAC  channels 
since  it  is  difficult  to  imagine  processes  that  would  make  one 
IRAC  channel  vary  without  the  other.  However,  this  method 
will  not  find  variables  in  cases  where  one  IRAC  channel 
is  compromised,  e.g.,  due  to  instrumental  effects,  and  the 
other  is  not. 

If  there  is  correlated  noise  between  the  two  channels  used 
for  the  Stetson  index,  especially  at  the  faint  end,  one  would 
expect  the  Stetson  index  to  be  correlated  with  source  brightness, 
as  seen  in,  e.g.,  Plavchan  et  al.  (2008b)  or  CHS01.  Figure  22 
shows  the  distribution  of  all  calculated  Stetson  indices  against 
mean  [3.6]  (for  the  standard  statistical  sample,  e.g.,  all  objects 
with  light  curves,  over  all  clusters,  fast  cadence  only).  Unlike 
the  analogous  figures  found  in,  e.g.,  Plavchan  et  al.  (2008b) 
or  CHS01,  here  we  have  no  substantial  change  in  the  bulk 
of  the  distribution  of  the  Stetson  index  toward  fainter  [3.6] 
magnitudes,  so  we  do  not  appear  to  have  correlated  noise 
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Figure  22.  Stetson  index  as  a  function  of  mean  [3.6]  in  magnitudes  for  all  clusters  for  the  standard  set  for  statistics  (YSOVAR-classic  fast  cadence  data).  The  lower 
panel  is  an  expanded  view  of  the  top  panel,  with  an  additional  gray  line  at  Stetson  index  =  0.  There  is  no  indication  of  correlated  noise  here,  as  would  be  the 
interpretation  of  a  change  in  the  distribution  of  Stetson  index  for  non- variables  as  a  function  of  brightness.  There  is  an  increase  in  frequency  of  large  Stetson  index 
values  for  brighter  [3.6],  but  that  is  likely  because  brighter  objects  are  more  likely  to  be  legitimately  young  cluster  members,  and  as  such  more  likely  to  be  variable 
(with  variability  correlated  between  the  two  IRAC  channels).  Note  that  essentially  no  variables  are  identified  fainter  than  [3.6]  ~  16,  consistent  with  our  observation 
that  those  light  curves  are  particularly  noisy. 


between  the  two  channels.  There  is  an  increase  in  frequency 
of  large  Stetson  index  values  for  brighter  [3.6],  but  that  is  likely 
because  brighter  objects  are  more  likely  to  be  legitimately  young 
cluster  members,  and,  as  such,  are  more  likely  to  be  variable 
(with  variability  correlated  between  the  two  IRAC  channels). 
Essentially  no  large  values  of  the  Stetson  index  are  identified  for 
objects  fainter  than  [3.6]  ~  16;  the  intrinsic  error  on  each  point  is 
sufficiently  large  that  these  objects  do  not  have  a  large  Stetson 
index.  (This  is  consistent  with  the  discussion  in  Section  4.2 
regarding  the  number  of  objects  in  our  data  set  falling  off  rapidly 
fainter  than  [3 .6]— [4.5]  ~  16.) 

The  specific  location  of  the  cutoff  between  variable  and  non¬ 
variable  can  be  unique  to  each  data  set,  as  it  is  affected  by  the 
sampling  length  and  rate  of  the  light  curves.  This  is  the  primary 
reason  behind  our  decision  to  calculate  statistics  over  only  the 
fast  cadence  window.  We  now  discuss  how  we  chose  this  cutoff 
value  for  the  Stetson  index.  The  left  panel  of  Figure  23  shows  a 
histogram  of  the  Stetson  indices  for  all  objects  in  the  standard  set 
for  statistics  having  at  least  five  points  in  both  1 1  and  12.  The  bulk 
of  the  distribution  about  0  are  the  non-variables,  and  that  part  of 
the  distribution  can  be  reasonably  well-fit  by  a  Gaussian.  There 
are  substantial  deviations  from  Gaussianity  toward  higher  values 
of  the  Stetson  index,  as  expected  for  a  population  of  identified 
variable  stars.  There  is  a  change  in  the  distribution  of  the  Stetson 
index  above  and  below  0.9;  from  where  the  distribution  deviates 
from  a  Gaussian  to  a  Stetson  index  of  0.9,  the  slope  in  the  middle 
panel  of  Figure  23  is  1.0,  and  from  a  Stetson  index  of  0.9  to 


~3,  the  slope  is  0.3.  Based  on  this,  we  take  0.9  as  the  cutoff 
for  variability  in  our  data  set.  The  value  of  0.9  corresponds  to 
about  6<t  for  the  Gaussian  fit  to  the  distribution. 

We  note  that  in  the  analogous  histograms  for  each  individual 
cluster,  each  typically  has  a  small  gap  in  the  Stetson  index 
distribution  at  ~0.9.  Orion,  however,  does  not,  and  Orion 
contributes  about  half  of  the  ~  11,000  viable  two-band  light 
curves  for  which  the  Stetson  index  appears  in  Figure  23. 

To  check  whether  the  Stetson  index  cutoff  of  0.9  is  sensible, 
we  conducted  a  series  of  Monte  Carlo  tests.  For  random 
light  curves  (with  a  Gaussian  distribution  of  points)  using 
the  same  time  sampling  as  the  real  data,  the  distribution 
of  Stetson  indices  is  well-described  by  a  Gaussian  with  a 
width  typically  of  0. 1-0.2,  comparable  to  the  left-hand  side 
of  Figure  23.  For  Orion,  where  we  have  a  reasonably  well- 
defined  set  of  members  and  non-members  (from  MCI  1;  not 
just  disked  and  non-disked,  but  confirmed  membership  lists), 
we  can  compare  the  distributions  of  Stetson  indices  for  the 
members  and  non-members.  In  Orion,  the  distribution  for 
non-members  is  generally  fairly  well-described  by  a  Gaussian 
centered  on  0,  but  there  is  a  small  “shoulder”  asymmetry  toward 
the  larger  Stetson  index  (likely  legitimate  field  variables  or  as 
of  yet  unidentified  members).  The  distribution  for  members,  in 
contrast,  is  not  well-described  by  a  Gaussian.  It  is  asymmetric 
with  a  substantial  excess  of  objects  with  high  Stetson  values. 
This  is  as  expected,  since  members  are  more  likely  to  have  large 
amplitude,  correlated  variability. 
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Figure  23.  Histograms  of  the  Stetson  indices  calculated  for  the  standard  set  for  statistics,  the  subset  of  which  has  sufficient  points  in  both  IRAC  bands.  Left:  linear 
histogram.  The  red  dotted  line  indicates  a  Gaussian  fit  to  the  histogram,  showing  deviations  from  Gaussianity  toward  higher  values  of  the  Stetson  index,  as  expected 
for  a  population  of  identified  variable  stars.  Middle:  a  zoomed  in  view  of  the  Stetson  values  between  0  and  5,  with  log  ordinate.  The  red  dotted  line  is  the  Gaussian 
fit  from  the  left  panel.  The  green,  dashed  lines  with  two  different  slopes  show  that  there  is  a  break  in  the  distribution  defining  our  cutoff  between  variable  (^0.9)  and 
non-variable  ( <0.9),  with  the  black  vertical  line  at  0.9.  A  value  of  0.9  corresponds  to  about  6a  for  the  Gaussian  fit  to  the  distribution.  Right:  histograms  of  the  Stetson 
indices  for  objects  in  our  standard  set  of  members  (blue  dashed  line)  and  likely  non-members  (black  solid  line).  While  this  division  is  imperfect,  the  distributions  cross 
at  ~0.9,  suggesting  that  relative  populations  of  members/non-members  is  the  dominant  effect  in  the  break  in  the  slope  of  the  entire  distribution  at  ~0.9.  We  conclude 
that  the  Stetson  index  cutoff  of  0.9  is  indeed  a  sensible  boundary  for  demarcating  variables  from  the  general  population. 

(A  color  version  of  this  figure  is  available  in  the  online  journal.) 


Similarly,  we  can  examine  the  distribution  of  the  Stetson 
index  for  our  standard  set  of  members  that  are  also  in  the 
standard  set  for  statistics  (and  having  sufficient  points  in  both 
bands)  and  compare  it  to  the  Stetson  index  distribution  for 
the  remaining  objects  not  selected  as  members  (but  still  in 
the  standard  set  for  statistics,  and  having  sufficient  points  in 
both  bands).  We  obtain  a  similar  result  in  the  right  panel  of 
Figure  23;  the  Stetson  distribution  for  members  crosses  that 
for  non-members  at  about  0.9,  or  perhaps  a  little  below  that 
level,  suggesting  that  this  division  is  the  dominant  cause  of 
the  break  in  the  entire  distribution  at  about  that  level.  Even  if 
our  separation  between  members  and  non-members  is  imperfect 
(which  it  certainly  is),  these  distributions  are  consistent  with  our 
selection  of  0.9  as  the  cutoff.  We  conclude  that  a  Stetson  index 
cutoff  of  0.9  is  a  sensible  boundary  for  our  data  set. 

While  the  objects  with  Stetson  indices  of  <0.4  have  a  very 
low  chance  of  having  legitimate  correlated  variability,  and  ob¬ 
jects  with  Stetson  indices  >0.9  have  a  high  likelihood  of  corre¬ 
lated  variability,  there  is  a  continuum  between  these  values. 
Objects  with  Stetson  indices  >0.4  and  <0.9  have  low- 
confidence  for  correlated  variability.  In  MCI  1,  we  took  a  differ¬ 
ent  Stetson  index  of  0.55  as  the  cutoff,  based  on  the  distribution 
of  Stetson  indices  for  that  particular  data  set.  As  such,  some 
of  the  identified  low-confidence  variables  may  have  changed 
between  this  and  the  initial  analysis.  We  also  note  that  the 
cutoff  in  Stetson  index  for  CSI  2264  is  very  different  (Cody 
et  al.  2014),  but  that  program  has  a  substantially  different  ob¬ 
serving  cadence  than  the  YSOVAR-classic  data  discussed  here. 
In  general,  the  appropriate  Stetson  index  cutoff  must  be  de¬ 
termined  for  the  individual  data  set,  and  there  is  no  universal 
value. 


5.2.  Chi-squared  Test 

A  second  method  to  identify  variables  is  a  chi-squared  test 
(x2),  which,  for  a  given  band,  is  given  by 

,  1  A  (mag,  -  mag)2 

— Vf — •  (3) 

»  =  1  1 

where  cr,  is  the  estimated  photometric  uncertainty  (corrected  as 
per  our  discussion  of  the  YSOVAR  noise  floor  in  Section  2.5). 

This  test  is  used  to  identify  objects  with  uncorrelated  vari¬ 
ability,  or  variability  in  only  one  band  (perhaps  because  data 
exist  in  only  one  band).  This  makes  the  x1  test  more  suscep¬ 
tible  to  instrumental  issues  affecting  only  one  band.  However, 
to  demonstrate  that  it  generally  does  a  good  job  of  recovering 
variables  with  large  Stetson  indices.  Figure  24  shows  the  distri¬ 
butions  of  xjt  and  /2,  as  a  function  of  Stetson  index.  For  those 
objects  in  the  standard  set  for  statistics  where  it  is  possible  to  cal¬ 
culate  both  x2  and  the  Stetson  index,  the  values  are  reasonably 
well  correlated  for  the  unambiguously  variable  objects.  Using 
this  plot,  we  find  that  a  limit  of  x2\  or  Xn  ~  5  is  an  appropriate 
conservative  cutoff  for  potential  variability  in  those  cases  where 
only  one  x2  can  be  calculated  (e.g.,  where  monitoring  in  only 
one  band  is  available). 

For  our  largest  data  set  (Orion),  there  are  thousands  of  light 
curves  that  meet  the  requirement  imposed  by  the  Stetson  index 
of  having  data  in  both  IRAC  channels.  However,  in  the  1 1 
smaller-field  YSOVAR-classic  data  sets,  imposing  such  a  two- 
band  restriction  typically  means  that  more  than  half  the  viable 
light  curves  would  be  discarded.  Thus,  the  x1  test  is  particularly 
useful  in  these  cases  where  only  one  band  is  available. 
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Figure  24.  Distributions  of  log  Xp  6|  (left)  and  log  X[4  5]  (right)  as  a  function  of  the  log  of  the  Stetson  index.  For  the  objects  within  the  standard  sample  for  statistics 
where  it  is  possible  to  calculate  both  x2  and  the  Stetson  index,  the  values  are  reasonably  well  correlated  for  the  variables.  The  vertical  line  is  at  a  Stetson  index  of  0.9, 
our  cutoff  for  selecting  variable  objects  based  on  the  Stetson  index.  On  the  basis  of  this  plot,  we  set  a  limit  of  X|2  6|  or  X|2  5.  ~  5  (horizontal  line)  as  the  cutoff  for 
potential  variability  in  those  cases  where  only  one  x2  can  be  calculated  (e.g.,  where  monitoring  in  only  one  band  is  available). 


Figure  25  shows  histograms  of  xj\  and  xji  f°r  the  standard 
sample  for  statistics,  as  well  as  for  the  subset  of  objects  for  which 
a  Stetson  index  can  be  calculated  (e.g.,  the  sample  used  in  the 
prior  figure),  and  the  smaller  subsample  of  objects  identified 
as  variable  using  the  Stetson  index  (Stetson  index  >0.9).  We 
fit  a  Gaussian  to  the  sharply  peaked  distributions,  and  found 
a  3cr  value  of  x2  <  4.5  in  all  cases.  The  bulk  of  the  x2 
distribution  for  which  the  Stetson  index  is  >0.9  also  has  x2  >  5. 
To  be  conservative,  we  thus  set  a  limit  of  xj\  °r  Xn  =  5 
for  identifying  a  candidate  variable  object.  In  the  ideal  case 
where  there  are  multiple  bands  of  monitoring  data,  one  could 
assess  each  light  curve  with  a  large  x2  but  small  Stetson  index 
for  physical  plausibility.  To  attempt  to  avoid  identifying  false 
variability  due  to  instrumental  effects,  the  light  curves  for  each 
of  these  potentially  variable  objects  will  be  examined  by  hand  in 
the  context  of  the  individual  cluster  analysis  to  come  in  separate 
papers. 


5.3.  Identifying  Periodic  Variables 

Finally,  as  in  MC 1 1 ,  there  are  still  legitimately  variable 
sources  within  our  standard  statistical  sample  that  fail  both 
the  Stetson  index  test  (perhaps  because  only  one  band  is 
available)  and  the  x2  test  (perhaps  because  the  amplitude  of 
variability  is  small  and  our  limits  for  identifying  variability 
were  conservative  by  design).  There  are  many  mathematical 
tools  available  for  identifying  periodic  behavior  in  an  unevenly 
sampled  time  series.  The  last  test  for  variability  we  run  here  is 
a  periodogram  analysis  using  the  NASA  Exoplanet  Archive 


Periodogram  Service34  (Akeson  et  al.  2013).  This  service 
provides  period  calculations  using  Lomb-Scargle  (LS;  Scargle 
1 982),  Box-fitting  Least  Squares  (BLS;  Kovacs  et  al.  2002),  and 
Plavchan  (Plavchan  et  al.  2008b)  algorithms.  These  methods  are 
varyingly  more  or  less  sensitive  to  periodic  behavior  shaped  like 
sinusoids  or  fiat-bottomed  transits,  and/or  may  be  less  sensitive 
to  periodic  behavior  appearing  in  addition  to  other  behavior, 
such  as  a  period  superimposed  on  a  long-term  trend.  The 
expected  periodic  variability  in  our  sample  includes  anything 
repeated,  from  a  sinusoidal-like  signal  originating  from  hot  or 
cool  spots  on  a  photosphere,  to  signals  characteristic  of  close 
binaries,  to  repeated  dips  in  the  signal  (like  “dippers”  or  AA  Tau; 
see,  e.g.,  MC11),  or  even  pulsations  (e.g.,  Morales-Calderon 
et  al.  2009). 

Specifically,  because  of  the  variety  of  expected  light  curve 
shapes,  and  the  weaknesses  inherent  in  any  of  these  methods  for 
finding  periodicity,  and  noting  the  approach  used  by  (and  results 
from)  McQuillan  et  al.  (2013a,  2013b),  we  also  calculated  the 
ACF  for  each  light  curve  as  a  check  on  repeated  patterns.  We 
linearly  interpolated  the  light  curve  onto  evenly  spaced  times, 
and  then  calculated  the  ACF  using  the  following  expression 
where  L  is  a  lag  in  days,  and  x  is  the  light  curve  (with 
elements  xf): 


ACFt(L)  =  ACFV(— L)  = 


Y!k=a  'fa  ~  x)(xk+L  -  x) 

k=0  (**  -  x> 


(4) 


34  http://exoplanetarchive.ipac.caltech.edu/cgi-bin/Periodogram/nph- 
simpleupload 


32 


The  Astronomical  Journal,  148:92  (46pp),  2014  November 


Rebull  et  al. 


I°9(X![3.6|)  l°9(X![4.5]) 


Figure  25.  Distributions  of  log  X|2  6|  (left)  ami  log  X|4  5|  (right)  for  all  objects  with  x2  values  (black  solid  line  histogram),  for  those  objects  with  x2  and  Stetson  indices 
(gray  solid  line  histogram),  and  for  those  objects  with  x2  and  a  Stetson  index  >0.9  (light  gray  filled  histogram).  (All  objects  shown  here  are  from  the  standard  set  for 
statistics.)  The  red  dotted  line  is  a  Gaussian  fit  to  the  corresponding  histogram.  The  3<r  values  corresponding  to  that  Gaussian,  converted  to  linear  /2,  are  indicated. 
We  identify  a  limit  of  Xp.g]  or  X[4  5|  ~  5  as  a  conservative  cutoff  for  potential  variability  in  those  cases  where  only  one  x2  can  be  calculated  (e.g.,  where  monitoring 
in  only  one  band  is  available).  That  limit  is  plotted  as  the  vertical  dashed  line.  We  take  objects  with  X[2  6|  or  X|2  5)  >  5  as  legitimately  variable. 

(A  color  version  of  this  figure  is  available  in  the  online  journal.) 


We  experimented  with  several  different  timescales  as  ob¬ 
tained  from  the  ACF,  and  settled  on  the  location  of  the  first 
peak,  providing  that  the  peak  was  above  an  ACF  value  of  0.2. 
For  those  objects  with  significant  periods,  this  coherence  time 
should  be  well-matched  to  the  period. 

We  looked  for  periods  in  light  curves  where  we  had  at  least  20 
points  (more  restrictive  than  our  standard  statistical  sample);  we 
ran  all  four  methods  (LS,  BLS,  Plavchan,  and  ACF)  on  not  just 
the  3.6  jint  and  4.5  gm  light  curves,  but  also,  where  possible, 
the  [3.6]— [4.5]  light  curves.  In  some  cases,  a  long-term  trend 
(astrophysical,  not  instrumental)  is  present  in  the  individual  II 
or  12  light  curve,  masking  a  periodic  signal,  but  the  color  exhibits 
the  periodic  signal.  We  looked  for  periods  between  only  0.05 
and  15  days,  given  the  overall  sampling  of  our  data,  and  we 
require  at  least  two  complete  periods  over  the  typically  ~40  day 
window  of  our  observations.  We  investigated  phased  light  curves 
for  those  periods  calculated  using  all  of  these  methods.  Based 
on  these  many  thousands  of  results,  we  concluded  that  LS  is  the 
best,  for  our  data  set,  for  finding  reliable,  plausible  periods.  BLS 
and  the  Plavchan  algorithm,  while  they  look  for  a  wider  variety 
of  shapes  of  signals,  struggle  with  light  curves  that  typically 
have  less  than  100  points,  as  ours  do;  the  ACF  approach  finds 
only  the  strongest  signals.  Typically,  if  the  LS  algorithm  found  a 
reliable  period,  those  other  three  approaches  found  comparable 
periods. 

Thus,  we  filtered  first  on  the  LS  results.  We  excluded  can¬ 
didate  periodic  objects  if  the  calculated  false  alarm  probability 
(FAP;  see,  e.g.,  Scargle  1982)  was  >0.03,  or  if  the  recovered 
period  was  >14.5  days  and  the  FAP  for  that  period  was  >0.01, 
or  if  the  period  was  <0. 1  day  (slightly  larger  than  the  lower 
limit  over  which  we  searched),  or  if  the  calculated  period  was 
exactly  15.0  days  (by  inspection  of  the  light  curves,  input  and 
phased,  a  returned  period  exactly  equal  to  the  upper  limit  of 


our  search  window  was  usually  indicative  not  of  a  true  period¬ 
icity,  but  instead  of  a  long-term  trend  in  the  data).  For  each  of 
the  remaining  objects,  we  investigated  the  phased  light  curve. 
Perhaps  unsurprisingly,  given  our  overall  FAP  cutoff  of  0.03, 
about  3%  of  the  surviving  candidate  periodic  light  curves  did 
not  produce  physically  plausible  phased  signals.  Those  objects 
were  omitted  from  the  final  set  of  periods,  and  will  be  identi¬ 
fied  as  such  in  the  corresponding  cluster  papers.  The  planned 
individual  cluster  papers  may  include  a  few  additional  periodic 
objects  not  automatically  identified  due  to  the  presence  of  outly¬ 
ing  photometry  points  which  mask  periodicity  unless  removed 
by  hand. 

We  proceed  to  include  all  of  the  periodic  objects  identified 
via  the  LS  algorithm  (excluding  the  candidate  periodic  objects 
as  described)  in  the  set  of  variable  objects,  even  if  they  fail  the 
other  (Stetson,  /2)  variability  tests.  Individual  objects  will  be 
discussed  in  the  papers  dedicated  to  each  individual  cluster,  but 
anywhere  from  1  to  15  objects,  typically  ^5,  were  added  to 
the  list  of  likely  variables  for  each  of  the  smaller-field  clusters. 
For  any  given  object,  we  wish  to  assign  a  single  period  to  that 
object.  We  take  any  period  derived  from  the  [3.6]  data  first  ([3.6] 
is  less  noisy  than  [4.5]),  then,  only  if  there  is  no  [3.6]  period 
of  sufficient  power,  we  take  the  period  derived  from  [4.5],  and 
finally,  if  no  other  period  of  sufficient  quality  is  available,  then 
we  take  that  derived  from  [3.6]-[4.5], 

A  preliminary  list  of  periodic  objects  in  Orion  appeared  in 
MC 1 1 .  The  approach  we  are  now  using  to  search  for  periods  is 
more  stringent  than  that  in  MC  1 1 .  Our  current  approach  recovers 
the  bulk  of  the  objects  that  MC11  reported  as  periodic,  but 
does  not,  for  example,  recover  the  objects  reported  as  having 
P  >  15  days.  A  complete  list  of  the  ~800  periodic  variables  in 
Orion  as  derived  from  this  YSOVAR  data  reduction  will  appear 
in  a  later  paper. 
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Figure  26.  Histograms  of  the  time  steps  between  observations  for  each  cluster.  In  all  cases,  the  time  between  observations  is  non-uniform  to  avoid  aliasing  for  a 
specific  period.  AFGL  490,  Mon  R2,  and  NGC  2264  have  a  lower  overall  sampling  rates  than  all  other  clusters. 


5.4.  Detection  Limits  for  Variability 

Above,  we  described  how  variable  sources  are  identified 
using  the  Stetson  test,  the  x 2  test,  and  a  search  for  periodic 
variability.  The  cut-off  values  for  those  tests  were  chosen  to 
yield  a  conservative  list  of  variable  sources  and  to  reject  those 
sources  where  the  variability  stems  mostly  from  instrumental 
artifacts.  We  now  present  Monte-Carlo  simulations  to  quantify 
how  much  variability  is  required  to  meet  those  criteria. 

Several  different  physical  effects  may  contribute  to  the 
observed  variability  and  this  can  lead  to  very  complex  patterns  in 
the  light  curve.  In  the  absence  of  a  theoretical  model  to  explain 
the  different  contributions,  we  concentrate  on  simple  analytical 
prescriptions  for  light  curves  so  as  to  gain  a  sense  of  the 
sensitivity  of  our  statistical  tools.  First,  we  consider  a  source  that 
has  two  states,  a  bright  and  a  faint  state.  The  light  curve  switches 
randomly  between  those  two  states.  We  simulate  light  curves  for 
a  different  fraction  of  time  spent  in  the  upper  state  (0%,  10%, 
20%,  30%,  40%,  and  50%)  and  we  expect  that  variability  is 
more  easily  found  for  sources  with  a  larger  amplitude  between 
the  two  states  and  an  equal  chance  to  find  the  source  in  each 


state.  Second,  we  simulate  sinusoidal  light  curves  with  different 
periods — 0.1-2  days  (in  steps  of  0.1  day)  and  2-20  days  (in 
steps  of  2  days).  We  add  Gaussian  noise  to  each  light  curve  and 
vary  the  ratio  of  the  signal  amplitude  and  the  noise  (from  0  to  1 0 
in  steps  of  0. 1 ,  where  a  relative  amplitude  of  0  means  a  constant 
light  curve  with  noise  only).  The  sensitivity  of  the  Stetson  test, 
the  x1  test,  and  the  Lomb-Scargle  periodogram  is  independent 
of  the  magnitude  of  a  source  (that  is  not  noise-dominated);  only 
the  relative  amplitude  of  the  signal  and  the  noise  plays  a  role. 
Thus,  it  is  equally  possible  to  detect  strong  variability  in  a  weak 
source  (with  large  photometric  uncertainties)  as  weak  variability 
in  a  bright  source  (with  small  photometric  uncertainties). 

For  each  grid  point  in  relative  amplitude  and  fraction  of 
time  in  the  upper  level  or  period,  for  each  band,  we  simulate 
10,000  light  curves  for  each  grid  point  for  each  band.  The  light 
curves  are  sampled  at  the  time  intervals  of  the  actual  YSOVAR 
observations.  Figure  26  shows  the  histogram  of  the  time  steps 
between  observations  for  each  cluster  (the  typical  min  and  max 
A /  were  given  above  in  Table  3).  For  all  clusters,  the  sampling 
is  non-uniform  to  avoid  aliasing  for  a  specific  period.  However, 
the  histograms  fall  in  two  groups,  with  AFGL  490,  Mon  R2,  and 
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Figure  27.  Representation  of  the  relative  time  steps  for  the  fast  cadence  observations,  in  a  similar  format  to  Figure  14,  with  a  “|”  denoting  a  time  step  (from  the  red 
sections  of  Figure  14).  AFGL  490,  Mon  R2,  and  NGC  2264  can  also  be  seen  here  as  having  a  lower  overall  sampling  rates. 
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Figure  28.  Efficiency  for  variability  detection.  The  left  panel  shows  simulations  for  light  curves  with  two  distinct  states  (bright  and  faint);  the  other  two  panels  show 
results  for  a  sinusoidal  light  curve.  The  contours  indicate  the  fraction  of  sources  in  that  part  of  the  diagram  that  return  x2  >  5.  In  the  first  two  panels,  colors  indicate 
the  probability  that  the  Stetson  index  >0.9  (and  note  that  the  color  scale  is  non-linear).  In  the  right  panel,  colors  indicate  the  simulated  period  that  was  recovered  in  a 
Lomb-Scargle  periodogram.  The  time  sampling  is  from  the  observations  of  L1688. 

(A  color  version  of  this  figure  is  available  in  the  online  journal.) 


NGC  2264  having  a  lower  overall  sampling  rate  than  the  other 
clusters.  Figure  27  is  another  representation  of  the  sampling  rate. 
In  the  style  of  Figure  14,  it  visually  represents  the  sampling  rates 
for  the  fast  cadence  monitoring.  Here,  AFGL  490,  Mon  R2,  and 
NGC  2264  can  also  be  seen  to  have  a  lower  overall  sampling 
rate  than  the  other  clusters.  We  ran  all  Monte  Carlo  simulations 
using  the  actual  time  sampling  from  LI 688,  one  of  the  clusters 
with  a  high  sampling  rate  and  Mon  R2,  a  representative  cluster 
with  a  lower  sampling  rate. 

Figure  28  shows  results  from  the  Monte  Carlo  simulations 
using  the  sampling  of  the  L1688  cluster.  The  left  panel  presents 
the  detection  efficiency  for  light  curves  from  an  object  with 
two  distinct  luminosities.  If  the  star  is  found  in  each  state  half 
the  time,  the  Stetson  test  will  identify  it  as  variable  in  almost 
all  cases  (99%)  if  the  step  size  is  at  least  three  times  larger 
than  the  noise  level.  Since  the  x2  test  uses  data  from  one  band 
only,  the  step  size  must  be  larger  (five  times  the  noise  level)  to 
reach  the  same  detection  efficiency.  If  the  star  spends  less  than 
20%  of  the  time  in  either  state,  there  is  a  reasonable  chance 
that  the  variability  will  not  be  found,  even  for  larger  step  sizes, 
since  the  sampling  might  find  only  few  data  points  in  this  state. 


For  a  periodic  light  curve  (middle  panel),  variability  is  again 
found  more  easily  in  the  Stetson  test  than  in  the  x2  test,  and 
the  period  of  the  variability  does  not  influence  the  detection 
efficiency,  since  those  two  tests  do  not  consider  the  time  ordering 
of  the  observed  data.  Using  the  LS  approach,  we  are  sensitive  to 
periods  ( P )  between  about  1  and  15  days.  Even  periods  where 
the  amplitude  (a)  of  the  light  curve  a  sin(27r  /  Pt)  is  only  twice  as 
large  as  the  noise  level  are  easily  detected,  almost  independent 
of  the  period  in  the  range  1-15  days.  Such  weak  signals  would 
not  necessarily  show  up  as  variable  in  the  Stetson  or  x2  test  (see 
the  middle  panel).  In  general,  the  region  where  the  tests  detect 
variability  in  some  light  curves  but  not  in  others  is  fairly  narrow. 

Figure  29  shows  the  same  plots  as  Figure  28  but  for  the  time 
sampling  of  Mon  R2,  one  of  the  clusters  with  a  lower  cadence 
in  the  observations.  The  general  shape  of  the  regions  in  the 
parameter  space  where  variability  can  be  detected  is  the  same, 
but,  due  to  the  lower  number  of  observations  and  the  larger  time 
span  between  observations,  a  larger  amplitude  is  required,  and 
we  are  less  sensitive  to  periods  below  about  two  days. 

To  explore  more  complicated  light  curves,  we  combined 
several  effects,  e.g.,  a  sinusoidal  periodicity  overlaid  on  a 
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Figure  29.  Same  as  in  Figure  28,  but  using  the  time  sampling  of  the  Mon  R2  observations.  The  different  time  sampling  effectively  means  that  slightly  increased 
signal-to-noise  ratios  are  required. 

(A  color  version  of  this  figure  is  available  in  the  online  journal.) 


long-term  trend.  In  these  cases,  the  strongest  effect  generally 
determines  how  the  variability  will  be  seen.  If  the  magnitude  of 
the  trend  is  large  and  the  amplitude  of  the  sine  wave  is  small,  then 
the  light  curve  will  be  marked  as  variable,  but  the  periodicity 
might  not  be  detected. 

Given  the  possible  complexity  of  real  light  curves,  it  is  not 
possible  to  cover  the  entire  parameter  space  with  Monte  Carlo 
simulations.  However,  the  scenarios  presented  here  show  that,  in 
general,  we  can  detect  periodicity  with  an  amplitude  just  twice 
the  level  of  the  noise,  variability  at  3-5  times  the  level  of  the 
noise  with  the  Stetson  test  and  about  6-10  times  the  level  of 
the  noise  for  single-band  light  curves  with  the  x 2  test.  Similar 
results  for  the  relative  sensitivity  of  the  Stetson  and  the  /2  test 
were  also  found  by  Flaherty  et  al.  (2013),  although  they  used 
different  cut  off  levels  than  this  work. 

As  seen  above  in  Figure  1 6  and  in  MC 1 1 ,  there  is  a  reasonably 
strong  correlation  between  the  mean  magnitude  of  a  source  and 
its  mean  error.  This  overall  relation  affects  our  ability  to  find 
variability  or  periodicity.  For  sources  brighter  than  13th  mag, 
the  total  uncertainty  is  dominated  by  the  error  floor  introduced  in 
Section  2.5.  For  those  brighter  sources,  we  can  detect  periodicity 
if  the  amplitude  is  larger  than  about  0.02  mag  and  variability 
if  it  is  larger  than  about  0.03-0. 1  mag  (with  the  exact  number 
depending  on  the  signal  shape). 

5.5.  Identifying  Variables:  Cryo-to-post-cryo 
(Six  to  Seven  Years) 

In  addition  to  the  variability  probed  on  the  YSOVAR  mon¬ 
itoring  timescale  of  ~40  days,  we  also  are  interested  in  the 
evidence  for  longer-timescale  variations  between  observations 
of  these  same  clusters  in  the  Spitzer  cryogenic  epoch  and  the 
post-cryo  (YSOVAR)  epochs.  To  identify  the  variables  in  this 
case,  for  every  object  with  a  light  curve,  we  can  compare  the 
average  measurement  from  the  earliest  cryo  era  (Table  4),  and 
the  mean  for  that  object  over  the  YSOVAR  standard  statistical 
sample  (Section  3.2). 

The  process  we  used  to  identify  the  long-term  variables  in 
each  cluster  is  shown  for  AFGL  490  in  Figure  30.  We  plot 
the  difference  between  the  cryo  and  post-cryo  measurements 
for  each  object  as  a  function  of  the  cryo  value  and  determine, 
for  each  cluster,  the  brightness  at  which  photometric  noise 
clearly  dominates.  This  faintness  limit  is  close  to  16th  mag, 
consistent  with  what  we  noted  in  Figures  19-22.  We  select 
this  limit  separately  for  each  cluster  and  consider  only  objects 


brighter  than  this  limit.  We  fit  a  Gaussian  to  the  histogram  of 
the  difference  between  the  cryo  and  post-cryo  measurements, 
allowing  the  zero  point  as  well  as  the  height  and  width  of 
the  Gaussian  to  be  free  parameters.  We  classify  as  long-term 
variables  all  objects  with  cryo  to  post-cryo  offset  further  than 
3er  from  the  peak  of  the  (fitted)  distribution.  All  objects  from 
the  standard  statistical  sample,  not  just  the  members  (or  only 
the  variables),  go  into  this  process  of  defining  the  width  of 
the  distribution.  The  fraction  of  objects  in  each  field  that  are 
classified  as  variable,  and  the  subset  of  variables  that  are  also 
cluster  members,  are  both  identified  after  the  long-term  variables 
are  identified.  A  summary  of  the  important  parameters  in  these 
analysis  steps  to  search  for  variables  over  this  long-term  baseline 
is  in  Table  7;  the  values  we  used  for  At,  the  time  lapse  between 
the  cryo  and  post-cryo  observations,  are  included  in  Table  4. 

6.  DISCUSSION 

In  this  section,  we  present  an  analysis  of  the  distribution 
of  rotation  rates  as  a  function  of  IR  excess,  evidence  (or  lack 
thereof)  for  transient  IR  excesses,  evidence  for  skews  over 
time  toward  more  brightening  or  fading  sources,  and  how  the 
long-term  variability  fraction  varies  as  a  function  of  cluster 
parameterization  (from  Section  4.1)  or  length  of  time  baseline 
sampled. 

6.1.  Periodic  Variables 

Our  YSOVAR  map  in  Orion  is  far  larger  than  the  other  cluster 
maps,  which  focus  on  the  most  embedded  (possibly  youngest) 
objects  in  these  clusters.  For  these  embedded  objects,  it  is  not 
generally  possible  to  obtain  a  rotation  period  from  ground-based 
optical  or  NIR  data  due  to  extinction.  Moreover,  for  stars  with 
more  significant  disks,  it  is  less  likely  that  the  IR  light  curve 
will  be  strictly  periodic.  Many  distinct  processes  can  contribute 
to  a  YSO’s  mid-IR  variability,  and  it  often  results  in  a  stochastic 
light  curve  (Cody  et  al.  2014). 

For  the  objects  for  which  we  can  derive  a  period,  we  would 
like  to  compare  our  values  to  those  from  the  literature  as  a  check 
on  our  methodology.  Since  our  clusters  are,  for  the  most  part, 
very  embedded,  there  are  not  very  many  known  periods  in  the 
literature.  As  discussed  in  Section  2. 1 .2,  there  are  many  periods 
for  Orion  and  NGC  2264,  typically  obtained  in  the  optical,  but 
with  some  values  from  the  NIR.  Parks  et  al.  (2014)  reports  on 
NIR  periods  from  objects  in  our  region  of  L 1 688;  there  are  other 
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Figure  30.  Illustration  of  the  process  used  to  identify  long-term  variables  in  the  clusters.  Top:  plot  of  the  difference  between  the  cryogenic-era  measurement  and  the 
mean  measurement  from  the  standard  set  for  statistics  against  the  cryogenic-era  measurement  (left:  [3.6),  right:  [4.5]).  Only  objects  brighter  than  15.5  were  used  (in 
this  cluster)  for  the  next  step.  Middle:  histograms  of  the  difference  in  magnitudes  for  objects  brighter  than  15.5.  The  solid  line  is  the  data;  the  dotted  line  is  a  Gaussian 
fit  to  the  histogram.  Based  on  this  fit,  we  retained  objects  for  which  the  difference  in  magnitudes  is  greater  than  3rr  away  from  the  peak  of  the  distribution  as  likely 
variables.  Bottom:  zoomed  in  view  of  the  central  portions  of  the  prior  two  histograms.  A  summary  of  this  analysis  searching  for  variables  over  this  long-term  baseline 
is  in  Table  7. 

(A  color  version  of  this  figure  is  available  in  the  online  journal.) 


Table  7 

Statistics  on  Variable  Objects  on  the  Longest  Timescales 


Cluster” 

1 1  Faint  Cutoff 
(mag) 

11  Center  ±a  b 
(mag) 

11  Vars/LC.sL 
(mag) 

11  Memb.vars/LCsd 

12  Faint  Cutoff 
(mag) 

12  Center  ±rr  b 
(mag) 

12  Vars/LCsc 
(mag) 

12  Memb.vars/LCs1 

AFGL  490 

15.8 

0.001  ±0.04 

86/907  =  0.09 

43/120  =  0.36 

15.8 

0.017  ±0.05 

112/904  =  0.12 

58/155  =  0.37 

NGC  1333 

16.0 

0.033  ±  0.05 

43/289  =  0.15 

24/115  =  0.21 

16.0 

0.035  ±  0.05 

82/393  =  0.21 

31/112  =  0.28 

Orion 

16.0 

0.008  ±  0.04 

738/5313  =0.14 

439/2326  =  0.19 

16.0 

0.013  ±0.04 

920/5921  =  0.16 

476/2426  =  0.20 

Mon  R2 

16.0 

0.020  ±  0.06 

57/431  =0.13 

32/170  =  0.19 

16.0 

0.031  ±0.11 

24/211  =0.11 

10/90  =  0.11 

GGD  12-15 

16.0 

0.014  ±0.05 

56/336  =  0.17 

35/127  =  0.28 

16.0 

0.046  ±  0.06 

44/400  =  0.11 

28/138  =  0.20 

NGC  2264 

15.8 

0.011  ±0.06 

47/365  =  0.13 

29/142  =  0.20 

15.8 

0.017  ±  0.09 

58/492  =  0.12 

27/199  =  0.14 

L1688 

15.5 

0.017  ±0.05 

34/264  =  0.13 

18/48  =  0.38 

15.5 

0.039  ±  0.06 

47/312  =  0.15 

13/43  =  0.30 

Serpens  Main 

15.5 

0.008  ±  0.04 

130/1112  =  0.12 

25/61  =0.41 

15.5 

0.038  ±  0.05 

135/1161  =0.12 

24/83  =  0.29 

Serpens  South 

15.5 

0.013  ±0.04 

59/587  =  0.10 

28/68  =  0.41 

15.5 

0.031  ±0.05 

80/706  =  0.11 

24/78  =  0.3 1 

IRAS  20050+2720 

15.8 

0.004  ±  0.05 

123/1372  =  0.09 

31/114=0.27 

15.8 

0.021  ±0.06 

112/1415  =  0.08 

32/106  =  0.30 

1C  1396A 

15.5 

-0.007  ±  0.03 

117/1122  =  0.10 

21/73=  0.29 

15.5 

0.003  ±  0.04 

185/2116  =  0.09 

44/175  =  0.25 

Ceph  C 

15.8 

0.004  ±  0.04 

51/531  =0.10 

27/71  =0.38 

15.8 

0.030  ±  0.05 

61/553  =  0.11 

33/82  =  0.40 

Notes. 

“  The  values  we  used  for  At,  the  time  laps  between  the  cryo  and  post-cryo  observations,  are  included  in  Table  4. 

b  Mean  II  or  12  cryo-to-postcryo  offset,  e.g.,  center  (peak  location)  and  cr  of  Gaussian  fit  to  distribution  of  differences  in  magnitudes  brighter  than  the  faint  cutoff. 
c  Number  of  identified  long-term  variables  over  all  of  the  objects  with  light  curves/number  of  light  curves  (for  all  objects)  =  long  term  variability  fraction  for  all  objects. 
d  Number  of  identified  long-term  variables  over  only  the  standard  set  of  members/number  of  light  curves  (for  members)  =  long  term  variability  fraction  for  members. 


37 


The  Astronomical  Journal,  148:92  (46pp),  2014  November 


Rebull  et  al. 


All  but  Orion  Just  Orion 


Figure  31.  YSOVAR-determined  period  in  days  against  the  literature  period  in 
days  for  (left)  NGC  2264,  LI 688,  and  IC  1396A,  and  (right)  Orion  alone.  Solid 
lines  indicate  a  1:1  match,  a  2:1  match,  and  a  1:2  match. 


literature  values  for  rotation  periods  of  objects  elsewhere  in 
LI 688,  beyond  our  monitored  region.  We  can  roughly  compare 
to  the  MIR  timescales  reported  in  Morales-Calderon  et  al.  (2009) 
forIC  1396A. 

Of  the  objects  with  periods  in  the  literature,  there  are 
~200  that  also  have  periods  derived  here  in  Orion  (exclud¬ 
ing  those  from  MCI  1,  since  those  were  derived  from  the  same 
observations  we  use  here),  and  an  additional  ~15  from  « 

NGC  2264,  L1688,  and  IC  1396A.  About  75%  of  those  have  ^ 

period  measurements  that  match  to  better  than  10%,  so  we  have  i= 
confidence  that  our  period-finding  approach  is  at  least  well- 
matched  to  those  in  the  literature.  Figure  3 1  plots  the  YSOVAR- 
determined  period  against  the  literature  period  for  those  objects 
where  it  is  possible.  The  clusters  that  are  not  Orion  (NGC  2264, 

LI 688,  and  IC  1396A)  are  plotted  separately  simply  because 
Orion  dominates  the  statistics,  and  it  is  useful  to  see  if  there  are 
good  matches  outside  Orion  as  well  as  within  Orion.  Three  of 
the  four  objects  from  the  three  smaller-field  clusters  that  are  not 
well-matched  to  the  YSOVAR-determined  period  are  close  to 
likely  harmonics,  and  the  periods  are  of  comparatively  low  qual¬ 
ity.  The  one  that  is  most  discrepant  is  from  NGC  2264,  SSTYSV 
064101.40+093408.1,  and  is  being  compared  to  a  period  from 
Lamm  et  al.  (2004).  Our  phased  light  curve  looks  correct  (for 
our  wavelength  and  epoch  of  observation).  Of  the  ~50  Orion 
periods  that  do  not  agree  to  10%  (out  of  ~200  Orion  period 
comparisons  total),  it  is  predominantly  the  case  (by  a  ratio  of 
3  to  1)  that  the  YSOVAR  period  is  longer  than  the  literature 
period.  About  60%  of  these  have  [3.6]— [8]  >  0.8.  All  but  five  of 
those  Orion  periods  were  optically  determined.  Because  we  are 
working  in  longer  wavelengths,  it  is  possible  that,  particularly 
in  those  cases,  we  may  be  sampling  a  different  location  in  the 
star-disk  system,  e.g.,  further  away  from  the  photosphere,  where 
Keplerian  rotation  periods  are  longer. 

A  relation  has  already  been  found  between  IR  excess  and 
rotation  rate  for  young  stars,  suggesting  that  IR  excess  and 
rotation  rate  are  related;  out  of  our  12  clusters,  this  relation 
has  been  found  in  Orion  (Rebull  et  al.  2006)  and  NGC  2264 
(Cieza  &  Baliber  2007).  In  the  1 1  smaller-field  clusters  (i.e., 
all  but  Orion),  there  are  ~350  stars  with  periods  measured 
from  YSOVAR  light  curves,  but  only  ~250  of  those  also  have 
cryogenic  Spitzer  measurements  at  3.6  and  8  /tin  from  which 
we  can  get  a  clear  indication  of  the  IR  excess  in  these  systems. 

There  are  ~800  stars  in  Orion  with  measured  YSOVAR  periods, 
but  only  ~430  have  [3.6]  and  [8]  measurements. 


All  but  Orion  All  incl.  Orion 
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Figure  32.  IR  excess  ([3.6]— [8])  vs.  log  (period  in  days)  for  objects  in  the 
YSOVAR  clusters,  using  the  periods  derived  from  the  YSOVAR  data  as 
described  in  the  text.  Left:  the  1 1  clusters,  excluding  Orion;  there  are  '-250 
objects.  Right:  all  12  YSOVAR  clusters,  including  Orion;  there  are  ~430  Orion 
objects  plus  the  ~250  objects  from  left  panel.  The  plots  are  similar,  both  to  each 
other  and  to  that  obtained  for  Orion  by  Rebull  et  al.  (2006),  despite  the  fact  that 
most  of  these  stars  are  on  average  thought  to  be  younger  than  those  in  Orion. 

All  but  Orion  All  incl.  Orion 


[3.6]-[8]  [3.6]-[8] 

Figure  33.  Cumulative  distributions  of  IR  excess  ([3.6]-|8J)  for  objects  with 
log(/5)  ^  0.25  (solid  line),  0.25  <  log(/J)  ^  0.75  (dotted  line),  and  0.75  <  log(R) 
(dashed  line).  Left:  the  1 1  smaller-field  clusters,  excluding  Orion.  Right:  all  12 
YSOVAR  clusters,  including  Orion.  The  distributions  are  significantly  different 
according  to  a  K-S  test. 

Figure  32  shows  the  relationship  between  IR  excess  (specif¬ 
ically  [3.6]— [8])  and  YSOVAR-derived  IR  periodicity  for  these 
sources.  In  both  cases,  there  is  a  gap  near  [3.6]— [8]  ~  0.8,  which 
divides  the  disk  candidates  (above  that  cutoff)  from  the  non-disk 
candidates  (below).  There  is  also  different  behavior  to  the  left 
and  right  of  log(R)  ~  0.25,  or  P  ~  1.8  days — excesses  do  not 
necessarily  imply  longer  periods,  but  a  star  with  a  longer  period 
is  more  likely  than  those  with  shorter  periods  to  have  an  IR  ex¬ 
cess.  Figure  33  shows  the  cumulative  distributions  of  [3 .6]— [8] 
for  the  same  two  panels  as  in  Figure  32,  for  three  different  bins 
of  log(R),  divided  at  log(R)  =  0.25  and  0.75  (1.78  days  and 
5.62  days,  respectively).  According  to  Kolmogorov-Smirnov 
(K-S)  tests,  the  distributions  of  [3.6]— [8]  are  significantly  dif¬ 
ferent  within  each  panel.  The  two  distributions  that  are  the  most 
similar  are  the  full  (all  12  clusters)  distributions  for  0.25  < 
log(R)  ^  0.75,  and  0.75  <  log(P);  the  probability  that  those 
populations  were  drawn  from  the  same  distribution  is  4%.  The 
probability  that  the  populations  were  drawn  from  the  same  dis¬ 
tribution  for  log(F)  ^  0.25  and  0.75  <  log(P)  (again,  for  all 
12  clusters)  is  ~  1 0  ~ 13 ;  for  log(R)  ^  0.25  and  0.25  <  log(P)  ^ 
0.75,  the  probability  that  the  populations  were  drawn  from  the 
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Figure  34.  Cumulative  distributions  of  log(P)  for  objects  with  [3.6j-[8J  > 
0.8  (solid  line),  and  [3.6]-[8]  ^  0.8  (dotted  line).  Left:  the  1 1  smaller-field 
clusters,  excluding  Orion.  Right:  all  12  YSOVAR  clusters,  including  Orion.  The 
distributions  are  significantly  different  according  to  a  K-S  test. 


same  distribution  is  ~10  l0.  Similarly,  Figure  34  shows  the  cu¬ 
mulative  distributions  of  log(P)  for  the  same  two  panels  as  in 
Figure  32,  for  two  different  bins  of  [3.6]-[8],  divided  at  [3.6]- 
[8]  =  0.8.  Again,  according  to  K-S  tests,  the  distributions  of 
log(P)  are  significantly  different  within  each  panel;  the  proba¬ 
bility  that  either  of  the  populations  were  drawn  from  the  same 
distribution  is  <10-17. 

The  plots  in  Figure  32  are  very  similar  to  that  obtained  for 
Orion  by  Rebull  et  al.  (2006),  and  that  by  other  investigators 
in  other  clusters,  despite  the  fact  that  optically  determined 
periods  were  used  there.  The  periods  derived  from  our  IR 
YSOVAR  data  may  be  photospheric  rotation  rates,  pulsation 
rates  (see,  e.g.,  Morales-Calderon  et  al.  2009),  or  inner  disk 
rotation  rates  (see,  e.g.,  Artemenko  et  al.  2012);  on  the  other 
hand,  for  those  clusters  where  there  are  periods  available,  we 
match  the  literature  reasonably  well,  and  the  literature  for  the 
most  part  is  using  optical  data  to  obtain  periods.  Therefore,  it 
seems  that  we  are,  in  most  cases,  not  sampling  much  different 
locations  in  the  star-disk  system.  However,  an  exception  could 
be  that  the  optical  and  IR  observations  are  sampling  two 
separate  places  whose  movements  are  locked  together,  such  as 
starspots  on  the  photosphere  and  stellar-magnetosphere-driven 
disk  disturbances  at  the  corotation  radius. 

It  is  also  surprising  that  the  results  for  the  aggregate  set  of 
clusters  are  so  similar  to  that  for  Orion,  because  the  clusters 
should  be  for  the  most  part  substantially  younger  than  Orion. 
This  could  imply  that  disk  locking  may  be  in  effect  at  even 
these  young  ages,  or  that  accretion-powered  stellar  winds  are 
the  dominant  mechanism  to  slow  these  objects.  However,  it 
is  likely  that  we  can  obtain  viable  periods  more  easily  for 
relatively  unobscured  stars,  e.g.,  with  these  periods,  we  are  also 
sampling  the  older  end  of  the  young  star  distributions  in  these 
clusters.  Little  is  known  about  many  of  the  objects  outside  of 
Orion  shown  here;  additional  study  of  the  individual  objects  will 
help  clarify  matters.  Individual  objects  will  be  discussed  in  the 
corresponding  YSOVAR  cluster  paper. 

6.2.  Disks  Do  Not  Vanish  or  Appear 

There  have  been  recent  reports  of  debris  disks  undergoing 
significant  short-term  changes;  both  Meng  et  al.  (2012)  and 
Melis  et  al.  (2012)  report  on  systems  that  change  significantly 
at  wavelengths  >10  /cm  over  timescales  of  years.  Our  objects 
are  much  younger  (a  few  million  years  rather  than  a  few  tens  or 


hundreds  of  million  of  years)  and  our  monitoring  wavelengths 
are  considerably  shorter.  However,  Rice  et  al.  (2012)  also  note 
that  nine  (36%)  of  the  stars  in  their  sample  of  young  stars  in 
Cygnus  OB7  (comparable  in  age  to  our  sample)  have  a  transient 
NIR  ( JHK )  excess.  We  can  use  this  first  look  at  our  data  to 
constrain  the  degree  to  which  IR  excesses  in  our  sample  vanish 
(or  appear)  on  the  timescales  of  years,  namely  between  the  cryo- 
era  observation  and  that  of  our  post-cryo  observations. 

Irrespective  of  whether  or  not  objects  have  been  identified 
as  variable  above,  we  compared  the  cryo-era  [3.6]— [4.5]  color 
with  the  maximum  and  minimum  [3.6]— [4.5]  color  obtained 
during  our  YSOVAR  (fast-cadence)  monitoring  (standard  set 
for  statistics,  the  subset  of  which  have  measurements  in  both 
channels).  Out  of  ~  11,000  objects  (cluster  members  as  well 
as  background  objects  included)  for  which  we  have  [3 .6]— [4.5] 
color  light  curves,  there  are  at  most  15  objects  that  seem  to 
have  legitimate  substantial  changes  to  the  [3.6]-[4.5]  color 
(changes  of  a  size  that  might  be  consistent  with  big  changes 
to  a  disk),  and  these  are  all  relatively  faint  ([3.6]  >  12  mag) 
objects.  At  most,  two  of  those  cases  have  an  IR  excess  that 
appears  to  be  possibly  transient  on  these  timescales,  so  at  most 
<0.02%  frequency  of  occurrence.  For  the  remaining  13  objects 
with  plausibly  real  changes  in  color,  the  disk  is  still  clearly 
present,  but  the  brightness  and  color  have  changed  substantially. 
Individual  objects  will  be  discussed  in  the  cluster  papers. 

6.3.  Brightening  as  Likely  as  Fading 

In  the  literature  (e.g.,  Giannini  et  al.  2009;  Antoniucci  et  al. 
2014),  constraints  have  been  placed  on  the  timescales  for 
brightening  or  fading  by  comparing  how  many  sources  are 
found  to  be  getting  brighter  or  getting  fainter  (e.g.,  for  each 
source,  given  the  two  epochs,  for  how  many  cases  is  the  second 
epoch  brighter  than  the  first,  and  for  how  many  cases  is  the  first 
epoch  brighter  than  the  second).  If  there  are  random  fluctuations 
in  brightness,  the  same  number  of  sources  should  get  brighter 
as  get  fainter.  If,  instead,  there  are  more  fading  sources  than 
brightening,  then  the  type  of  variability  may  be  characterized 
by  a  short  rise  and  a  long  fall. 

We  can  make  a  similar  comparison  among  our  long-term 
variable  sample.  To  reduce  scatter  from  noise,  we  consider  only 
the  standard  set  of  members  (Section  3.1),  and  consider  only 
those  objects  identified  as  long-term  variables  independently 
in  both  [3.6]  and  [4.5]  (Section  5.5).  Figure  35  compares  the 
numbers  of  these  remaining  sources  that  become  brighter  at  both 
[3.6]  and  [4.5]  with  those  that  become  fainter  at  both  [3.6]  and 
[4.5]  for  1 1  of  the  clusters;  errors  are  approximated  by  simple 
Poisson  counting  statistics.  A  significantly  larger  number  of 
variables  are  found  in  Orion,  with  essentially  equal  numbers  of 
sources  getting  brighter  as  getting  fainter.  For  the  smaller-field 
clusters,  the  numbers  of  sources  that  brighten  is  less  consistently 
equal  to  the  number  of  sources  that  fade,  but  there  are  also 
far  fewer  sources  to  count.  Summing  up  the  1 1  smaller-field 
clusters,  there  are  107  (±10)  sources  that  brighten  and  88  (±9) 
sources  that  fade;  these  numbers  are  within  2a  of  each  other, 
but  both  numbers  have  to  be  extended  ~  1  a  toward  each  other. 
Fitting  a  line  to  the  points  in  Figure  35,  including  Orion,  results 
in  a  line  of  slope  1.05  ±  0.12,  and  an  intercept  of  —2.8  ±  3.3. 
This  is  consistent  with  no  difference  between  the  brightening 
and  fading  sources,  but  with  only  a  ~1ct  possibility  that  there 
are  slightly  more  brightening  sources.  With  or  without  Orion, 
there  are  very  similar  numbers  brightening  and  fading. 

In  the  case  of  several  other  papers  in  the  literature,  they 
were  looking  for  much  more  significant  variability  than  we 
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Figure  35.  Plot  of  the  number  of  members  that  become  fainter  vs.  the  number 
that  become  brighter  in  both  [3.6J  and  [4.5]  between  the  cryo  epoch  and  the 
YSOVAR  epoch.  Error  bars  are  approximated  by  Poisson  statistics.  The  Orion 
point  is  to  the  far  upper  right,  at  ( 1 56,  153),  with  errors  of  ~  1 2  in  each  direction. 
The  gray  line  is  the  unity  relation.  There  are  similar  numbers  of  objects  that 
become  brighter  as  become  fainter;  see  the  text. 


are  finding  (e.g.,  FUors  and  EXors).  They  found  internally 
consistent  patterns  of  sharp  rises  and  long  falls  in  brightness. 
For  our  sample,  we  are  apparently  finding  more  varieties  of 
variability,  even  on  the  long-term,  such  that  the  timescales 
average  out  over  six  to  seven  years  to  have  no  significant  bias 
toward  brightenings  or  fadings. 

6.4.  Long-term  Variability  Fractions 

Flaving  identified  the  long-term  variables  above  (Section  5.5) 
for  each  cluster,  we  can  look  at  the  fraction  of  objects  that 
are  variable  on  the  longest  timescales  we  sample.  Figure  36 
shows  the  fraction  of  variables  for  all  objects  in  each  field 
(the  standard  set  for  statistics),  including  both  members  and 
likely  field  objects.  The  x  axis  is  our  parameterization  of  the 
relative  fractions  of  embedded  sources  (see  Section  4.1),  the 
ratio  of  Class  I/total  number  of  members.  Figure  36  thus  is 
somewhat  incongruous  in  that  the  standard  set  of  members  is 
used  in  the  x  axis,  but  the  y  axis  includes  everything  in  the  FOV. 
About  10%-20%  of  all  objects  with  light  curves  are  identified 
as  variable  in  the  long-term.  However,  as  can  be  seen  in  the 
figure,  the  variability  fraction  of  everything  in  the  field  is  not 
a  strong  function  of  the  Class  I/total  ratio — the  slopes  of  the 
best-fit  lines  in  Figure  36  are  both  small.  (The  slopes  of  these 
lines  are  given  in  the  figure  caption.)  Calculating  Pearson’s 
correlation  coefficient  also  suggests  that  there  are  no  significant 
correlations  shown  in  Figure  36.  We  expect  that  the  fraction  of 
stars  that  are  variable  should  be  a  function  of  Galactic  latitude 
(see  approximate  Galactic  latitudes  listed  in  Table  1 ),  because 
the  fraction  of  sources  that  are  background/foreground  stars 
will  be  higher  in  the  Galactic  plane,  so  the  fraction  of  cluster 
members  will  be  higher  out  of  the  Galactic  plane,  and  since 
any  young  star  (cluster  member)  is  more  likely  to  be  variable. 


there  should  be  a  higher  fraction  of  variable  objects  at  higher 
Galactic  latitude.  Indeed,  the  Pearson’s  correlation  coefficient 
suggests  that  the  fraction  of  long-term  variables  for  all  objects 
in  each  field  is  strongly  correlated  with  the  absolute  value  of  the 
Galactic  latitude. 

Figure  37  recasts  Figure  36  by  using  the  long-term  variability 
fraction  of  only  the  standard  set  of  members,  rather  than  every¬ 
thing  in  the  field.  This  time,  there  is  a  significant  correlation;  the 
higher  the  fraction  of  embedded  members,  the  higher  the  frac¬ 
tion  of  long-term  variables.  If  a  cluster  has  more  sources  that 
are  embedded  and  likely  to  be  actively  accreting  and  interacting 
with  their  circumstellar  material,  it  also  has  more  sources  that  are 
variable  on  timescales  of  years.  The  slopes  of  the  best-fit  lines 
shown  in  Figure  37  are  given  in  the  figure  caption.  Pearson’s 
correlation  coefficient  suggests  that  there  is  a  correlation  here, 
and  it  is  stronger  in  II  than  12  (consistent  with  the  calculated 
slopes). 

Serpens  South  has  the  highest  fraction  of  the  most  embedded 
sources;  to  test  if  it  is  providing  a  significant  “lever  arm”  on  this 
fit,  we  fit  the  remaining  points  omitting  Serpens  South,  which 
does  indeed  weaken  the  correlation,  as  can  be  seen  in  the  figure. 

We  tested  the  correlations  seen  in  Figures  36  and  37  using 
a  variety  of  other  SED  class-based  parameterizations,  and  we 
found  similar  results — there  is  no  significant  correlation  of  the 
long-term  variable  fraction  of  all  sources  in  the  field  with  any 
SED  class-based  parameterization,  and  there  is  a  correlation 
of  the  long-term  variable  fraction  of  the  members  with  any 
parameterization  chosen  that  uses  fractions  of  various  SED 
classes  (or  groups  of  classes).  This  correlation  between  the 
fraction  of  members  that  are  variable  on  the  longest  timescales 
and  the  parameterization  of  “embeddedness”  seems  robust.  We 
expected  that  young  stars  were  more  likely  to  be  variable  than 
field  stars.  Moreover,  we  have  found  that  within  the  category 
of  young  stars,  for  a  higher  fraction  of  embedded  sources  (a 
higher  fraction  of  presumably  younger  sources),  we  find  a  higher 
fraction  of  long-term  variables.  This  is  consistent  with  what  has 
been  found  in  individual  clusters  (e.g.,  NGC  2264  (Cody  et  al. 
2014)  and  L1688  (Gunther  et  al.  2014)). 

The  difference  between  fractions  of  embedded  objects  among 
these  clusters  is  not  the  only  potentially  significant  factor  for 
this  set  of  observations.  Here,  we  sample  timescales  of  ~4.5  to 
~7.5  yr.  If  the  amplitude  of  variability  of  the  members  increases 
as  the  time  baseline  increases,  we  expect  more  members  to  be 
selected  as  variable  in  the  long-term,  and  thus  a  higher  long¬ 
term  variability  fraction  as  the  time  baseline  increases.  However, 
Figure  38  shows  this  relationship  between  the  variability  frac¬ 
tion  and  the  time  between  epochs  of  observation  (from  Table  4). 
The  best-fit  lines  and  correlation  coefficients  are  consistent  with 
no  significant  effect  on  the  variability  fraction  as  a  function  of 
timescale  sampled.  Serpens  South,  because  it  was  observed  (in¬ 
deed,  discovered;  Gutermuth  et  al.  2008b)  comparatively  late 
in  the  Spitzer  mission,  samples  the  shortest  timescales.  If  the 
Serpens  South  point  is  omitted  from  the  fit,  the  slopes  become 
slightly  steeper,  but  there  is  still  no  significant  correlation  from 
the  correlation  coefficients. 

Scholz  (2012),  working  in  the  A"  band,  finds  that  the  longer 
one  monitors  a  cluster,  the  larger  the  amplitude  of  variability, 
on  timescales  of  years.  That  work  specifically  investigated  the 
amplitude  of  the  change  in  magnitude,  for  just  one  cluster  at  a 
time,  and  looked  at  the  change  of  the  range  of  that  distribution  of 
amplitudes  as  a  function  of  time  step.  For  the  times  we  sampled 
on  these  longest  timescales,  analogous  plots  do  not  show  a 
significant  change  in  the  median  or  the  top  quartile  or  the  90th 
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Figure  36.  Fraction  of  long-term  variables  from  the  standard  set  for  statistics,  as  a  function  of  the  ratio  of  Class  I/total  objects,  for  all  objects  in  the  field,  for  [3.6]  (top) 
and  [4.5]  (bottom),  for  each  cluster,  as  indicated.  The  range  in  the  y  axis  is  set  to  match  the  range  needed  for  the  next  figure.  Error  bars  are  calculated  assuming  Poisson 
counting  statistics;  Orion  has  by  far  the  most  sources,  and  so  has  by  far  the  smallest  error  bars.  The  gray  lines  are  the  best-fit  line,  using  errors  in  both  directions  for 
each  point.  The  long-term  variability  fraction  for  everything  in  the  field  is  relatively  constant  in  this  plot — the  slopes  of  these  lines  are,  for  [3.6],  —0.14  ±  0.04,  and 
for  [4.5],  —0.06  ±  0.06.  Correlation  coefficients  suggest  that  there  is  no  significant  correlation  here  ([3.6]:  r  —  — 0.41 ;  [4.5J:  r  =  —0. 10). 
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Figure  37.  Similar  to  Figure  36,  except  that  the  y  axis  is  the  fraction  of  long-term  variable  members.  The  slopes  of  the  gray  best-fit  solid  lines  are,  for  [3.6],  0.66  ± 
0.19,  and  for  [4.5],  0.45  ±0.18.  Correlation  coefficients  are  consistent  with  a  correlation  ([3.6]:r  =  0.58;  [4.5]:r  =  0.41).  Since  Serpens  South  has  the  largest  fraction 
of  embedded  objects,  we  also  tested  fitting  this  relation  omitting  this  point,  and  this  does  weaken  the  correlation.  The  dashed  lines  are  these  best  fit  values:  for  [3.6], 
0.42  ±  0.15,  and  for  [4.5],  0.31  ±  0.14.  (Correlation  coefficients  are  [3.6]:r  =  0.45;  [4.5]:r  =  0.42)  The  long-term  variability  fraction  increases  significantly  as  the 
degree  of  embeddedness  increases  (to  the  right). 


percentile.  We  tried  using  the  most  stringent  set  of  objects  (only 
those  in  the  standard  set  of  members,  standard  set  for  statistics, 
that  had  light  curves  in  both  IRAC  channels);  we  still  did  not  find 
this  effect.  For  this  analysis,  Scholz  (2012)  was  only  working  in 
p  Oph,  over  a  larger  area  than  we  were,  using  slightly  different 
wavelengths,  and  a  larger  range  of  timescales,  the  maximum  of 
which  (~2000  days)  is  comparable  to  the  minimum  At  (~4.5  yr 


~1600  days)  we  consider  here.  These  could  account  for  the 
observed  differences. 

7.  CONCLUSIONS 

We  present  in  this  paper  the  data  collection  and  reduction 
for  the  YSOVAR  programs,  representing  nearly  800  hr  of 
Spitzer  time  studying  the  variability  of  young  stars  in  12 
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Figure  38.  Long-term  variable  fraction  for  the  standard  set  of  members  as  a  function  of  the  time  difference  in  years  between  the  cryo-epoch  and  YSOVAR  fast-cadence 
epoch.  The  slopes  of  the  gray  best-fit  lines  are,  for  [3.6],  0.04  ±  0.02  (r  =  -0. 16),  and  for  [4.5],  0.02  ±  0.02  (r  =  -0.0004).  Excluding  Serpens  South  as  an  outlier, 
the  slopes  of  the  gray  dashed  best-fit  lines  are  0.07  ±  0.02  and  0.04  ±  0.02,  respectively.  Correlation  coefficients  for  these  options  are  consistent  with  there  being  no 
significant  correlation  in  either  case  ([3.6]:r  =  0.35;  [4.5]:r  =  0.23).  All  of  this  evidence  is  consistent  with  no  significant  effect  of  the  timescale  on  the  fraction  of 
long-term  variables. 


different  clusters  (AFGL  490,  NGC  1333,  Orion,  Mon  R2, 
GGD  12-15,  NGC  2264,  L1688,  Serpens  Main,  Serpens  South, 
IRAS  20050+2720,  IC  1396A,  and  Ceph  C).  We  also  describe 
the  assembly  of  broad  collections  of  ancillary  data  for  these 
clusters.  There  are  ~29,000  unique  objects  of  any  sort  matched 
to  39,000  [3.6]  or  [4.5]  light  curves  in  the  YSOVAR  data  set. 

The  goals  of  the  broader  YSOVAR  program  include  the 
following:  to  obtain  the  first  extensive  mid-infrared  time  series 
photometry  of  young  stars  to  help  reveal  the  structure  of  the  inner 
disk  region  of  YSOs,  provide  new  constraints  on  accretion  and 
extinction  variability,  assess  timescales  of  mid-IR  variability 
from  seconds  to  years,  identify  new  young  eclipsing  binaries, 
help  identify  new  very  low-mass  substellar  members  of  the 
surveyed  clusters,  constrain  the  short-  and  long-term  stability 
of  hot  spots  on  the  surfaces  of  YSOs,  and  determine  rotational 
periods  for  objects  too  embedded  for  such  monitoring  in 
the  optical. 

In  this  paper,  we  set  the  stage  for  several  planned  papers. 
We  establish  here  not  only  the  data  reduction  approach,  but 
also  define  the  standard  sample  on  which  we  calculate  statistics 
(the  fast  cadence  data,  where  there  are  at  least  five  points  per 
light  curve),  and  a  standard  sample  of  members  (the  union  of 
all  IR-selected  members  and  X-ray-selected  members),  with  a 
provision  for  adding  members  identified  in  other  ways  (literature 
or  variability  itself). 

We  use  three  mechanisms  to  identify  variables  in  the  standard 
set  for  statistics  (fast  cadence  data) — the  Stetson  index  (calcu¬ 
lated  using  both  IRAC  channels),  the  x  2  test  (calculated  for  each 
channel  individually),  and  searching  for  significant  periodicity 
(working  on  light  curves  with  at  least  20  points,  using  primarily 
the  LS  approach,  independently  on  [3.6],  [4.5],  and  [3.6]— [4.5]). 
Based  on  simulations,  for  these  YSOVAR  data,  we  find  that  we 
are  sensitive  to  a  broad  range  of  timescales  and  amplitudes.  If 
the  star  is  found  in  one  of  two  states  half  the  time,  the  Stetson 
test  will  identify  it  as  variable  in  almost  all  cases  (99%)  if  the 


step  size  between  states  is  at  least  three  times  larger  than  the 
noise  level.  If  the  star  spends  less  than  20%  of  the  time  in  either 
state,  then  there  is  a  reasonable  chance  that  the  variability  will 
not  be  found,  even  for  larger  step  sizes  between  states,  since 
the  sampling  might  catch  only  few  data  points  in  this  state.  In 
general,  we  can  detect  periodicity  with  an  amplitude  twice  the 
level  of  the  noise,  variability  at  3-5  times  the  level  of  the  noise 
with  the  Stetson  test  and  6-10  times  the  noise  for  single-band 
light  curves  with  the  x2  test. 

We  also  identified  variables  on  the  longest  timescales  possible 
using  our  data  set,  timescales  of  six  to  seven  years,  using  a  fourth 
method  of  identifying  variability.  By  comparing  measurements 
taken  early  in  the  Spitzer  mission  with  the  mean  from  our 
YSOVAR  campaign  (the  standard  set  for  statistics),  we  can 
identify  those  objects  that  have  significantly  changed  between 
then  and  now.  We  show  that  the  overall  fraction  of  everything  in 
each  field  that  varies  on  these  longest  timescales  is  independent 
of  the  ratio  of  Class  I/total  members  in  each  cluster.  However, 
the  fraction  of  members  in  each  cluster  that  are  variable  on 
these  longest  timescales  is  a  function  of  the  ratio  of  Class  1/ 
total  members  in  each  cluster,  such  that  clusters  that  have  a 
higher  fraction  of  Class  I  objects  also  have  a  higher  fraction 
of  long-term  variables.  We  find  no  dependence  of  the  fraction 
of  members  in  each  cluster  that  are  variable  on  these  longest 
timescales  with  the  time  step  between  the  observations  (between 
the  cryogenic  and  post-cryogenic  observations).  Among  the 
most  reliable  of  the  long-term  variables,  we  find  no  strong 
preference  for  brightening  or  fading  over  these  timescales. 

We  find  periods  from  our  data  in  ~  1 1 00  objects,  ~800  of 
which  are  in  Orion  alone.  About  650  of  those  have  data  in 
both  3.6  and  8  pm  (~430  of  which  are  in  Orion),  enabling  us  to 
compare  [3.6]— [8]  versus  log(P)  in  a  fashion  similar  to  that  found 
previously  in  Orion  (Rebull  et  al.  2006)  and  NGC  2264  (Cieza 
&  Baliber  2007).  Very  similar  results  are  obtained — excesses 
do  not  necessarily  imply  longer  periods,  but  a  star  with  a  longer 
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period  is  more  likely  than  those  with  shorter  periods  to  have  an 
IR  excess.  This  is  somewhat  surprising,  in  that  ( 1 )  the  periods  are 
determined  from  the  IR,  not  the  optical,  as  was  done  previously; 
and  (2)  the  clusters  besides  Orion  are  thought  to  be  substantially 
younger  than  Orion,  suggesting  that  disk  locking  may  be  in 
effect  at  even  these  young  ages.  However,  it  is  likely  that  we 
can  obtain  viable  periods  more  easily  for  relatively  unobscured 
stars,  e.g.,  the  older  end  of  the  young  star  distributions  in  these 
clusters.  Little  is  known  about  many  of  the  objects  outside  of 
Orion  shown  here;  additional  study  of  the  individual  objects  will 
help  clarify  matters. 

There  have  been  recent  reports  of  debris  disks  changing 
on  timescales  of  years;  both  Meng  et  al.  (2012)  and  Melis 
et  al.  (2012)  report  on  systems  that  change  significantly  at 
wavelengths  >  1 0  jim  over  timescales  of  years.  Our  objects  here 
are  younger  and  our  monitoring  wavelengths  are  considerably 
shorter.  Out  of  ~  11,000  objects  (cluster  members  as  well  as 
foreground/background  objects  included)  for  which  there  is 
an  essentially  simultaneous  YSOVAR  measurement  in  both 
[3.6]  and  [4.5],  there  are  at  most  15  objects  that  seem  to  have 
legitimate  substantial  changes  to  the  [3 .6]— [4.5]  color  (changes 
of  a  size  that  might  be  consistent  with  a  disk  appearing/ 
disappearing),  and  these  are  for  the  most  part  relatively  faint 
([3.6]  >  12  mag)  objects.  At  most,  two  of  those  cases  have 
an  IR  excess  that  appears  to  be  possibly  transient  on  these 
timescales,  so  at  most,  <0.02%  frequency  of  occurrence.  For 
the  remaining  13  objects  with  plausibly  real  changes  in  color, 
the  disk  is  still  clearly  present,  but  the  brightness  and  color  have 
changed  substantially. 

Details  of  individual  objects  of  interest  in  each  of  the  clusters 
will  appear  in  the  forthcoming  YSOVAR  papers. 
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APPENDIX  A 
NAMING  CONVENTION 

In  our  initial  data  release  (MC11),  we  used  a  naming  con¬ 
vention  following  the  IAU  naming  standards  using  an  acronym 
(Initial  Spitzer  Orion  YSOVAR:  ISOY)  followed  by  the  J2000 
coordinates. 

For  the  final  version  of  our  YSOVAR  catalog,  discussed 
here,  the  IAU-registered  acronym  is  SSTYSV,  for  Spitzer 
Space  Telescope,  YSOVAR.  We  again  follow  it  with  the  J2000 
coordinates.  Individual  objects  in  this  catalog  need  not  be 
confirmed  young  stars,  but  simply  have  a  light  curve  in  this  data 
set.  Detailed  data  tables  of  cluster  members  for  each  cluster 
will  be  presented  in  the  individual  cluster  papers,  and  it  is  our 
intention  to  deliver  a  final  catalog  of  every  object  with  a  light 
curve  to  IRSA  for  general  distribution. 

APPENDIX  B 

SED  CLASSES 

B.I.  Background 

In  the  context  of  understanding  the  evolution  of  young 
stars  from  a  very  embedded  state  to  a  less  embedded  state, 
the  community  has  chosen  to  parameterize  objects  based  on 
the  slope  of  the  SED  (see,  e.g.,  Wilking  et  al.  2001).  This 
classification  is  tied  to  the  empirical  shape  of  the  SED.  Very 
embedded  (presumably  very  young)  objects  will  have  SEDs  that 
peak  at  long  wavelengths;  as  the  object  sheds  its  natal  cocoon, 
the  peak  of  the  SED  moves  to  shorter  wavelengths.  Objects 
with  substantial  circumstellar  disks  will  emit  more  energy  in  the 
IR  (and  longer  wavelengths)  than  in  the  optical.  Objects  with 
less  substantial  disks  will  have  most  of  their  energy  emitted  in 
wavelengths  shorter  than  the  NIR,  but  there  will  be  additional 
energy  contributions  in  the  IR  (and  longer  wavelengths)  from 
the  circumstellar  dust  and/or  debris,  i.e.,  the  SED  has  an  IR 
excess.  The  most  embedded  phase  is  referred  to  as  Class  0,  then 
proceeding  (based  on  SED  shape)  through  Class  I  (rising  SED), 
Flat  (flat  SED),  Class  II  (falling  SED  with  an  IR  excess),  and 
finally  Class  III  SEDs,  which  are  objects  with  photospheric  or 
near-photospheric  SEDs,  with  little  or  no  IR  excess,  typically 
identified  as  young  via  other  means  such  as  X-rays  or  Ha 
emission.  Sometimes  additional  classes  such  as  transition  disks 
are  added  near  the  end  of  this  sequence.  Classical  T  Tauri 
stars  (CTTSs)  are  often  identified  with  Class  IIs,  and  weak- 
lined  T  Tauri  stars  (WTTSs)  are  often  identified  with  Class  Ills. 
Nomenclature  is  difficult  and  inconsistent  across  the  literature; 
see  Evans  et  al.  (2009a)  for  a  discussion  of  terminology. 
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We  need  to  establish  at  least  an  internally  consistent  definition 
of  the  SED  classes  so  that  we  can  investigate  trends  as  a  function 
of  SED  class  as  a  proxy  for  age.  The  reliability  of  the  translation 
between  SED  class  and  age  has  been  discussed  at  length  in  the 
literature  (and  will  continue  to  be  discussed  in  the  future);  other 
factors  such  as  inclination  and  multiplicity  may  play  a  large  role. 
In  the  context  of  our  work,  we  wish  to  establish  an  internally 
consistent  approach  that  can  be  calculated  for  all  sources  in  the 
YSOVAR  fields. 

B.2.  Definition 

In  order  to  assemble  our  SEDs  for  each  object,  we  include  all 
the  data  described  in  Section  2  between  U  and  25  jxm.  For  the 
SED  classes,  only  the  IR  bands  are  relevant. 

In  order  to  define  an  internally  consistent  placement  of  the 
YSOVAR  objects  into  SED  classes,  in  the  spirit  of  Wilking  et  al. 
(2001),  we  define  the  near-  to  mid-IR  (2-24  gm)  slope  of  the 
SED,  a  =  d  log  XFx/d  log  X,  where  a  >  0.3  for  a  Class  I,  0.3 
to  —0.3  for  a  flat-spectrum  source,  —0.3  to  —1.6  for  a  Class  II, 
and  <  —  1 .6  for  a  Class  III.  For  each  of  the  objects  in  our  sample, 
we  performed  a  simple  least  squares  linear  fit  to  all  available 
photometry  (only  detections,  not  including  upper  or  lower  limits, 
but  including  archival  and  literature  data)  as  observed  between 
2  and  24  pm,  inclusive.  We  included  the  mean  obtained  from 
the  standard  set  for  statistics  of  the  light  curve  in  the  SED,  in 
addition  to  the  cryo-era  Spitzer  points.  (This  means  that  there 
could  be  more  than  one  point  contributing  to  the  fit  at  3.6  p m, 
one  from  the  cryo  era  and  one  from  the  mean  of  the  YSOVAR 
light  curves.)  Formal  errors  (either  from  individual  single-epoch 
measurements  or  the  mean  and  standard  deviation  from  the 
IRAC  light  curves)  on  the  infrared  points  are  so  small  as  to  not 
affect  the  fitted  SED  slope.  However,  if  the  mean  is  calculated 
over  more  of  the  light  curve  than  our  standard  set  for  statistics 
(fast  cadence)  sample,  the  mean  may  be  different  enough,  in  the 
cases  of  sparse  SEDs,  that  the  class  may  change  to  an  adjacent 
class;  these  cases  will  be  noted  in  the  cluster  papers  where 
relevant.  The  linear  fit  is  performed  on  the  observed  SED,  e.g., 
no  reddening  corrections  are  applied  to  the  observed  photometry 
before  fitting. 

In  the  literature,  the  precise  definition  of  a  can  vary,  or 
the  distribution  of  slopes  and  classes  can  vary  (e.g.,  Sung 
et  al.  2009  lists  Class  I,  II,  II/III,  pre-transition  disk,  transition 
disk  categories  for  NGC  2264),  which  may  result  in  different 
classifications  for  certain  objects.  Classification  via  our  method 
is  provided  specifically  to  enable  comparison  within  this  paper 
(and  to  other  YSOVAR  papers)  via  internally  consistent  means. 
Our  classification,  since  it  is  based  on  observed  SED,  is  possible 
for  all  objects  (not  just  those  identified  as  YSO  candidates). 
The  formal  classification  definition  puts  no  lower  limit  on  the 
colors  of  Class  III  objects  (thereby  including  those  with  SEDs 
resembling  bare  stellar  photospheres,  and  allowing  for  other 
criteria  such  as  X-ray  brightness  to  define  youth).  The  SED 
slopes  and  classes  for  all  of  the  sources  of  interest  will  appear  in 
the  individual  cluster  papers.  Histograms  of  the  relative  fractions 
of  each  class  for  the  standard  set  of  members  for  each  cluster 
appear  in  Figure  17. 

B.3.  Including  or  Ignoring  the  24  pm  Point 

In  terms  of  aggregate  statistics  over  all  data  from  all  12 
clusters  (members  and  non-members  together),  we  can  constrain 
the  fraction  of  objects  that  change  class  depending  on  whether  or 
not  the  24  pm  point  is  included  in  the  fit.  There  are  about  21,000 


objects  with  light  curves  over  all  12  clusters  for  which  an  SED 
slope  can  be  fit  between  2  and  8  pm.  Out  of  those  ~2 1,000, 
there  are  only  about  1760  with  MIPS-24  detections  (not  limits), 
so  only  about  8%  of  the  sources  are  affected.  Admittedly,  those 
sources  that  are  detectable  at  24  pm  are  the  ones  with  rising 
SEDs  and  thus  are  statistically  more  likely  to  be  true  cluster 
members  than  a  source  selected  at  random  from  the  map.  Out 
of  the  --1760  with  MIPS-24  detections,  ~72%  of  them  do  not 
change  class  when  the  24  pm  point  is  included  in  the  SED  fit. 
The  class  bins  of  slopes  are  relatively  large  and  thus  relatively 
insensitive  even  to  a  point  at  the  far  red  end  of  the  SED.  Of  the 
~28%  that  do  change  class,  ~85%  move  only  one  step,  to  an 
adjacent  class.  As  expected,  there  is  a  bias,  when  including  the 
24  pm  point,  to  move  the  objects  to  a  more  positive  slope,  e.g., 
toward  the  more  embedded  end  of  the  sequence;  of  the  ones  that 
change  class  at  all,  ~61%  move  one  step  earlier  (toward  more 
embedded,  not  necessarily  in  age)  in  the  sequence,  ~24%  of 
those  move  one  step  later  (toward  less  embedded,  not  necessarily 
in  age)  in  the  sequence. 

We  conclude  that  it  does  not  make  a  significant  difference  for 
the  overwhelming  majority  of  the  sources  if  one  uses  the  slope 
between  2  and  8  pm  or  the  slope  between  2  and  24  pm.  (Note 
that  this  is  fitting  all  available  points  between  these  values,  not 
a  simple  comparison  of  the  two  end  points.)  For  any  sources 
of  interest  in  which  the  class  might  change  depending  on  the 
inclusion  of  the  24  pm  point,  they  will  be  noted  in  the  individual 
cluster  papers. 

B.4.  Comparison  to  G09 

G09  provides  placement  into  SED  classes  as  part  of  the  data 
presented  there,  and  the  same  algorithm  has  been  applied  to 
our  entire  cryogenic-era  catalog  (Section  2.7).  These  classes  are 
identified  based  on  dereddened  colors,  and  are  only  provided 
for  objects  thought  to  be  young  stars.  Objects  that  are  not 
thought  to  be  young  stars  are  identified  as  other  things  such 
as  “PAH  emission  source”;  objects  missing  bands  such  that  the 
classification  scheme  cannot  be  run  are  also  not  classified.  The 
SED  classes  we  are  using  here  are  based  on  fits  to  the  observed 
SED,  not  the  dereddened  SED,  and  are  obtained  for  any  source, 
regardless  of  its  true  underlying  nature.  To  constrain  the  degree 
to  which  we  might  be  introducing  a  bias  by  fitting  the  observed 
SED,  we  can  compare,  for  some  objects,  the  class  obtained  by 
G09  and  by  our  mechanism  above. 

It  is  important  to  note  that  the  classes  provided  by  G09  are 
different  than  the  classes  we  use  here.  G09  has  Class  0, 1,  II,  and 
II/III,  but  no  Flat  class.  Here,  we  do  not  have  Class  0;  we  have 
Class  I,  Flat,  Class  II,  and  Class  III.  Over  all  12  clusters,  there 
are  ~3500  objects  for  which  both  classifications  are  available; 
75%  of  those  do  not  change  class,  even  with  the  different  bins 
that  are  defined.  Out  of  the  ~3500,  about  14%  are  identified  as 
being  from  a  later  (less  embedded)  class  than  G09,  and  about 
1 1%  are  identified  as  being  from  an  earlier  (more  embedded) 
class  than  G09. 

G09  also  dereddens  the  SEDs  before  placing  them  in  classes; 
we  are  fitting  observed  SEDs.  G09  does  this  to  avoid  a  reddening 
bias  toward  youth — as  discussed  in  Muench  et  al.  (2007),  the 
Xj— [3.6]  color  can  be  affected  by  Ay  ~  40,  though  it  takes 
Ay  >  200  to  affect  the  [5.8]— [24]  color.  However,  in  our  case, 
because  we  are  fitting  all  available  points  between  Ks  and  [24], 
not  just  subtracting  two  points  in  the  SED,  the  influence  of 
reddening  on  this  overall  slope  in  the  best  case  (where  there 
are  Ks,  four  bands  of  IRAC,  and  one  MIPS  band),  simulations 
suggest  that  we  would  need  Aj  ~  11  or  A v  ~  40  before  a 


44 


The  Astronomical  Journal,  148:92  (46pp),  2014  November 


Rebull  et  al. 


Class  III  object  would  be  misclassified  as  a  Class  II.  Admittedly, 
this  is  a  best-case  scenario,  where  the  SED  is  well-populated. 
For  a  source  without  MIPS  24  but  just  Ks  through  [8],  we  find 
that  we  need  Aj  ~  6  or  Ay  ~  21  before  a  Class  III  object 
would  be  misclassified  as  a  Class  II. 

Our  objects  discussed  in  YSOVAR  have  to  be  bright  enough 
to  get  good  quality  light  curves  at  3.6  and  4.5  gm,  which  is 
effectively  brighter  than  16th  mag  (see  Figures  19  and  20),  which 
means  that  our  sources  cannot  generally  have  extremely  high 
extinction.  Out  of  the  objects  for  which  we  have  a  Gutermuth- 
derived  value  for  A%  and  a  light  curve,  ~6%  have  A  k  >2 
(which  corresponds  to  Ay  >  6),  and  ~1%  have  A*-  >3.5 
(which  corresponds  to  Ay  >  11).  We  conclude  that  any  bias 
toward  young  objects  is  likely  not  substantial  in  our  data  set. 

We  have  opted  to  use  most  often  our  SED  class  definition 
as  described  above  for  internal  comparison  and  consistency. 
However,  when  discussing  the  Class  II/Class  I  ratio  reported 
in  G09,  we  are  using  those  classes  from  G09.  In  Sections  2.1.2 
and  4.1,  we  discuss  the  Class  II/Class  I  ratio  derived  in  a  very 
similar  fashion  as  G09,  using  the  G09  classes,  but  only  for  those 
objects  with  light  curves  in  the  standard  set  for  statistics.  This 
enables  at  least  some  comparison  back  to  G09  and  related  works. 

APPENDIX  C 

OBSERVATION  FOOTPRINTS 

In  support  of  the  discussion  in  Section  2.3,  this  Appendix 
includes  the  approximate  sky  coverage  for  a  summed  up  image 
consisting  of  all  epochs  of  the  YSOVAR  observations  for  each 
of  the  clusters  (except  for  AFGL  490,  which  is  included  in 
the  main  body  of  the  text).  In  each  case,  footprint  outlines 
are  superimposed  on  a  reverse  grayscale  image  of  the  cluster 
at  4.5  /x m  obtained  during  the  cryogenic  mission.  The  thicker 
blue  solid  line  is  3.6  /xm  and  the  thicker  red  dashed  line  is 
4.5  /xm.  If  there  is  substantial  field  rotation  during  the  YSOVAR 
campaigns,  a  single  epoch  of  observation  is  also  indicated  by 
the  thinner  blue  solid  and  red  dashed  lines,  with  the  difference 
between  the  single  epoch  and  the  larger  polygon  due  to  (ecliptic 
latitude  dependent)  field  rotation  effects;  see  Section  2.3.  North 
is  up  and  east  is  to  the  left  in  each  case.  The  distance  between 
the  farthest  north  and  farthest  south  coverage  here  is  noted  in 
each  caption.  The  relevant  Chandra  coverage  in  each  cluster  is 
shown  as  a  yellow  polygon. 
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